Meta-analysis of avian responses and butterfly eyespots

1 Setting ups

1.1 Load packages and custum function

Load packages and custum function. Custom function is used for calculating effect size (lnRR) and its variance from raw data.

Code

pacman::p_load(ape,
               GoodmanKruskal,
               DT,
               dtplyr,
               ggokabeito,
               ggcorrplot,
               here, 
               knitr,
               tidyverse,
               patchwork,
               metafor,
               miWQS,
               multcomp,
               orchaRd,
               parallel)

# function for calculating effect size (lnRR)
effect_lnRR <- function(dt) {
  # replace 0 in "C_mean", "T_sd", "C_sd", "C_proportion" with each minimum values
  # if proportion too extreme value, replace minimum value (only "T_proportion")

  dt <- dt %>%
  mutate(
    C_mean = ifelse(C_mean == 0, 0.04, C_mean),
    T_sd = ifelse(T_sd == 0, 0.01, T_sd),
    C_sd = ifelse(C_sd == 0, 0.05, C_sd),
    C_proportion = ifelse(C_proportion == 0, 0.01, C_proportion),
    T_proportion = ifelse(T_proportion < 0.01, 0.01, T_proportion),
    T_proportion = ifelse(T_proportion == 1, 0.9, T_proportion)
  )
    
  # copy dataset for adding effect size and its variation (lnRR /lnRR_var) column
  dt1 <- dt %>%
    mutate(
      lnRR = NA,
      lnRR_var = NA
    )

  for (i in seq_len(nrow(dt1))) {
    Tn <- dt1$Tn[i]
    Cn <- dt1$Cn[i]
    T_mean <- dt1$T_mean[i]
    C_mean <- dt1$C_mean[i]
    T_proportion <- dt1$T_proportion[i]
    C_proportion <- dt1$C_proportion[i]
    T_sd <- dt1$T_sd[i]
    C_sd <- dt1$C_sd[i]
    Response <- dt1$Response[i]
    Measure <- dt1$Measure[i]
    Study_design <- dt1$Study_design[i]
    Direction <- dt1$Direction[i]

    # continuous data - using escalc() (metafor package)
    if (Response == "continuous" & Study_design == "independent") {

      effect <- escalc(
        measure = "ROM",
        n1i = Tn, n2i = Cn,
        m1i = T_mean, m2 = C_mean,
        sd1i = T_sd, sd2i = C_sd,
        var.names = c("lnRR", "lnRR_var")
      )

      dt1$lnRR[i] <- effect$lnRR
      dt1$lnRR_var[i] <- effect$lnRR_var
    } else if (Response == "continuous" & Study_design == "dependent") {

      effect <- escalc(
        measure = "ROMC",
        ni = (Tn + Cn) / 2,
        m1i = T_mean, m2 = C_mean,
        sd1i = T_sd, sd2i = C_sd,
        ri = 0.5,
        var.names = c("lnRR", "lnRR_var")
      )

      dt1$lnRR[i] <- effect$lnRR
      dt1$lnRR_var[i] <- effect$lnRR_var
    }

    # proportion data (no sd values)
    else if (Response == "proportion1" & Study_design == "independent") {
      
    T_proportion <- replace(
        T_proportion, Direction == "Decrease",
        (1 - T_proportion[Direction == "Decrease"])
      )
      C_proportion <- replace(
        C_proportion, Direction == "Decrease",
        (1 - C_proportion[Direction == "Decrease"])
      )

      # transform proportion value
      asin_trans <- function(proportion) {
        trans <- asin(sqrt(proportion))
        trans
      }

      T_proportion <- asin_trans(T_proportion)
      C_proportion <- asin_trans(C_proportion)

      # calculate lnRR and lnRR variance
      lnRR_pro1 <- log(T_proportion / C_proportion)
      lnRR_var_pro1 <- (1 / sqrt(8))^2 * (1 / (T_proportion^2 * Tn) +
        1 / (C_proportion^2 * Cn))

      dt1$lnRR[i] <- lnRR_pro1
      dt1$lnRR_var[i] <- lnRR_var_pro1

    } else if (Response == "proportion1" & Study_design == "dependent") {
      
      T_proportion <- replace(
        T_proportion, Direction == "Decrease",
        (1 - T_proportion[Direction == "Decrease"])
      )
      C_proportion <- replace(
        C_proportion, Direction == "Decrease",
        (1 - C_proportion[Direction == "Decrease"])
      )

      # transform proportion value
      asin_trans <- function(proportion) {
        trans <- asin(sqrt(proportion))
        trans
      }

      T_proportion <- asin_trans(T_proportion)
      C_proportion <- asin_trans(C_proportion)

      # calculate lnRR and lnRR variance
      lnRR_pro1 <- log(T_proportion / C_proportion)
      lnRR_var_pro1 <- (1 / sqrt(8))^2 * (1 / (T_proportion^2 * Tn) +
        1 / (C_proportion^2 * Cn)) -
        2 * 0.5 * sqrt((1 / sqrt(8))^2 * (1 / (T_proportion^2 * Tn))) *
          sqrt((1 / sqrt(8))^2 * (1 / (C_proportion^2 * Cn)))

      dt1$lnRR[i] <- lnRR_pro1
      dt1$lnRR_var[i] <- lnRR_var_pro1
    }

    # proportion data (exist sd values)
    else if (Response == "proportion2" & Study_design == "independent") {
      
      T_proportion <- replace(
        T_proportion, Direction == "Decrease",
        (1 - T_proportion[Direction == "Decrease"])
      )
      C_proportion <- replace(
        C_proportion, Direction == "Decrease",
        (1 - C_proportion[Direction == "Decrease"])
      )
      
      # transform proportion mean value
      asin_trans <- function(proportion) {
        trans <- asin(sqrt(proportion))
        trans
      }

      T_SD <- sqrt(T_sd^2 / (4 * (T_proportion) * (1 - (T_proportion))))
      C_SD <- sqrt(C_sd^2 / (4 * (C_proportion) * (1 - (C_proportion))))

      T_proportion <- asin_trans(T_proportion)
      C_proportion <- asin_trans(C_proportion)

      # calculate lnRR and lnRR variance
      lnRR_pro2 <- log(T_proportion / C_proportion)
      lnRR_var_pro2 <- (T_SD)^2 * (1 / (T_proportion^2 * Tn)) +
        (C_SD)^2 * (1 / (C_proportion^2 * Cn))

      dt1$lnRR[i] <- lnRR_pro2
      dt1$lnRR_var[i] <- lnRR_var_pro2

    } else if (Response == "proportion2" & Study_design == "dependent") {
      
      T_proportion <- replace(
        T_proportion, Direction == "Decrease",
        (1 - T_proportion[Direction == "Decrease"])
      )
      C_proportion <- replace(
        C_proportion, Direction == "Decrease",
        (1 - C_proportion[Direction == "Decrease"])
      )
      
      # transform proportion mean value
      asin_trans <- function(proportion) {
        trans <- asin(sqrt(proportion))
        trans
      }

      T_SD <- sqrt(T_sd^2 / (4 * (T_proportion) * (1 - (T_proportion))))
      C_SD <- sqrt(C_sd^2 / (4 * (C_proportion) * (1 - (C_proportion))))

      T_proportion <- asin_trans(T_proportion)
      C_proportion <- asin_trans(C_proportion)

      # calculate lnRR and lnRR variance
      lnRR_pro2 <- log(T_proportion / C_proportion)
      lnRR_var_pro2 <- (T_SD)^2 * (1 / (T_proportion^2 * Tn)) +
        (C_SD)^2 * (1 / (C_proportion^2 * Cn)) -
        2 * 0.5 * sqrt((T_SD)^2 * (1 / (T_proportion^2 * Tn))) *
          sqrt((C_SD)^2 * (1 / (C_proportion^2 * Cn)))

      dt1$lnRR[i] <- lnRR_pro2
      dt1$lnRR_var[i] <- lnRR_var_pro2
    }
  }

  return(dt1)
}

1.2 Prepare datasets for analysis

Read datasets and phylogeny

First, we used a predator dataset to check whether we need to consider phylogeney. If so, we will use predator and prey datasets for following meta-analysis.

Code

dat_all <-  read.csv(here("data/all_15082023.csv"), header = T, fileEncoding = "CP932")
dat_pred <- filter(dat_all, Dataset == "predator")
dat_prey <- filter(dat_all, Dataset == "prey")

datatable(dat_pred, caption = "Predator datasets", options = list(pageLength = 10, scrollX = TRUE))

Table 1 - Predator datasets

Code

trees <- read.nexus(here("data/bird_phy.nex"))
plot(trees[1])

Figure 1— Phylogenetic tree of bird species

Calculate lnRR and lnRR variance and prepare datasets for analysis

Code

# predator dataset
dt_pred <- effect_lnRR(dat_pred)
dt_pred$Obs_ID <- 1:nrow(dt_pred)
dt_pred$Log_diameter <- log(dt_pred$Diameter_pattern)
dt_pred$Log_area <- log(dt_pred$Area_pattern)
dt_pred$Log_background <- log(dt_pred$Area_background)

# prey dataset
dt_prey <- effect_lnRR(dat_prey)
dt_prey$Obs_ID <- 1:nrow(dt_prey)
dt_prey$Log_diameter <- log(dt_prey$Diameter_pattern)
dt_prey$Log_area <- log(dt_prey$Area_pattern)
dt_prey$Log_background <- log(dt_prey$Area_background)

# all dataset
dt_all <- effect_lnRR(dat_all)
dt_all$Obs_ID <- 1:nrow(dt_all)
dt_all$Log_diameter <- log(dt_all$Diameter_pattern)
dt_all$Log_area <- log(dt_all$Area_pattern)
dt_all$Log_background <- log(dt_all$Area_background)

hist(dt_all$lnRR)

Code

hist(dt_all$lnRR_var)

Code

# use vcalc to calculate variance-covariance matrix
VCV_pred <- vcalc(vi = lnRR_var, 
                  cluster = Cohort_ID, 
                  subgroup = Obs_ID, 
                  data = dt_pred)

VCV_prey <- vcalc(vi = lnRR_var, 
                  cluster = Cohort_ID, 
                  subgroup = Obs_ID, 
                  data = dt_prey)

VCV <- vcalc(vi = lnRR_var, 
             cluster = Cohort_ID,
             obs = Obs_ID,
             rho = 0.5,
             data = dt_all)

summary(dt_pred)

   Study_ID          Cohort_ID         Shared_control_ID       Year     
 Length:117         Length:117         Length:117         Min.   :1957  
 Class :character   Class :character   Class :character   1st Qu.:2009  
 Mode  :character   Mode  :character   Mode  :character   Median :2009  
                                                          Mean   :2005  
                                                          3rd Qu.:2013  
                                                          Max.   :2016  
                                                                        
   Country          Data_source        Data_location      Bird_species      
 Length:117         Length:117         Length:117         Length:117        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
 Bird_common_name   Prey_species       Prey_common_name     Response        
 Length:117         Length:117         Length:117         Length:117        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
   Measure           Direction         Detail_result      Note_result       
 Length:117         Length:117         Length:117         Length:117        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
 Treatment_stimulus Number_pattern Diameter_pattern  Area_pattern    
 Length:117         Min.   :1.00   Min.   : 2.20    Min.   :  3.801  
 Class :character   1st Qu.:1.00   1st Qu.: 4.50    1st Qu.: 15.904  
 Mode  :character   Median :2.00   Median :20.00    Median : 87.715  
                    Mean   :1.94   Mean   :15.29    Mean   :254.622  
                    3rd Qu.:2.00   3rd Qu.:25.00    3rd Qu.:489.697  
                    Max.   :5.00   Max.   :31.28    Max.   :490.874  
                                                                     
 Area_background   Area_ratio      Shape_pattern      Control_stimulus  
 Min.   :  225   Min.   :0.01383   Length:117         Length:117        
 1st Qu.: 1000   1st Qu.:0.03623   Class :character   Class :character  
 Median : 3200   Median :0.07434   Mode  :character   Mode  :character  
 Mean   : 5851   Mean   :0.09281                                        
 3rd Qu.:13175   3rd Qu.:0.11257                                        
 Max.   :26400   Max.   :0.33071                                        
                                                                        
 Note_stimulus       Type_prey          Shape_prey         Diet_bird        
 Length:117         Length:117         Length:117         Length:117        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
   Bird_sex           Bird_age               Tn              Cn       
 Length:117         Length:117         Min.   : 4.00   Min.   : 4.00  
 Class :character   Class :character   1st Qu.: 8.00   1st Qu.: 8.00  
 Mode  :character   Mode  :character   Median : 9.00   Median : 9.00  
                                       Mean   :14.73   Mean   :14.26  
                                       3rd Qu.:15.00   3rd Qu.:15.00  
                                       Max.   :89.00   Max.   :89.00  
                                                                      
     T_mean              C_mean             MT_mean         
 Min.   :   0.0004   Min.   :   0.0400   Min.   :-556.1000  
 1st Qu.:   0.7656   1st Qu.:   0.5079   1st Qu.:  -6.8415  
 Median :  16.1000   Median :   9.3000   Median :   0.7544  
 Mean   : 119.1540   Mean   : 110.1728   Mean   : -43.7981  
 3rd Qu.: 153.5000   3rd Qu.: 154.2000   3rd Qu.:   0.9522  
 Max.   :1180.4878   Max.   :2115.0000   Max.   :   0.9996  
 NA's   :16          NA's   :16          NA's   :81         
    MC_mean                T_sd               C_sd            T_proportion    
 Min.   :-2115.0000   Min.   :  0.0100   Min.   :   0.0500   Min.   :0.01000  
 1st Qu.:  -13.6000   1st Qu.:  0.5183   1st Qu.:   0.2859   1st Qu.:0.05258  
 Median :    0.5543   Median : 16.9814   Median :   4.7434   Median :0.27152  
 Mean   : -121.8786   Mean   : 82.2725   Mean   :  78.2657   Mean   :0.30879  
 3rd Qu.:    0.9174   3rd Qu.: 74.7486   3rd Qu.:  56.2141   3rd Qu.:0.49500  
 Max.   :    0.9574   Max.   :731.8482   Max.   :1052.1749   Max.   :0.82609  
 NA's   :81           NA's   :16         NA's   :16          NA's   :75       
  C_proportion     Study_design         Dataset          Mean_median       
 Min.   :0.04262   Length:117         Length:117         Length:117        
 1st Qu.:0.08262   Class :character   Class :character   Class :character  
 Median :0.26477   Mode  :character   Mode  :character   Mode  :character  
 Mean   :0.29854                                                           
 3rd Qu.:0.50000                                                           
 Max.   :0.79618                                                           
 NA's   :75                                                                
 T_mean_median       C_mean_median       Type_of_variance_statistic
 Min.   :   0.0004   Min.   :   0.0000   Length:117                
 1st Qu.:   0.7656   1st Qu.:   0.6559   Class :character          
 Median :  12.0000   Median :   9.3000   Mode  :character          
 Mean   : 104.2470   Mean   :  96.8310                             
 3rd Qu.: 104.0000   3rd Qu.:  61.0000                             
 Max.   :1180.4878   Max.   :2115.0000                             
                                                                   
      T_se                C_se          T_variance_q1  T_variance_q3
 Min.   :  0.00000   Min.   :  0.0000   Min.   :77.5   Min.   :246  
 1st Qu.:  0.09462   1st Qu.:  0.0903   1st Qu.:77.5   1st Qu.:246  
 Median :  4.88000   Median :  1.5000   Median :77.5   Median :246  
 Mean   : 25.83973   Mean   : 26.7360   Mean   :77.5   Mean   :246  
 3rd Qu.: 24.68328   3rd Qu.: 19.8747   3rd Qu.:77.5   3rd Qu.:246  
 Max.   :182.92683   Max.   :372.0000   Max.   :77.5   Max.   :246  
 NA's   :20          NA's   :20         NA's   :116    NA's   :116  
 C_variance_q1 C_variance_q3     Note                lnRR         
 Min.   :29    Min.   :36    Length:117         Min.   :-2.15769  
 1st Qu.:29    1st Qu.:36    Class :character   1st Qu.:-0.07696  
 Median :29    Median :36    Mode  :character   Median : 0.08779  
 Mean   :29    Mean   :36                       Mean   : 0.27498  
 3rd Qu.:29    3rd Qu.:36                       3rd Qu.: 0.54880  
 Max.   :29    Max.   :36                       Max.   : 4.16777  
 NA's   :116   NA's   :116                                        
    lnRR_var             Obs_ID     Log_diameter       Log_area    
 Min.   :0.0008295   Min.   :  1   Min.   :0.7885   Min.   :1.335  
 1st Qu.:0.0107871   1st Qu.: 30   1st Qu.:1.5041   1st Qu.:2.767  
 Median :0.0333636   Median : 59   Median :2.9957   Median :4.474  
 Mean   :0.1249595   Mean   : 59   Mean   :2.4180   Mean   :4.574  
 3rd Qu.:0.1655345   3rd Qu.: 88   3rd Qu.:3.2189   3rd Qu.:6.194  
 Max.   :1.3885984   Max.   :117   Max.   :3.4430   Max.   :6.196  
                                                                   
 Log_background  
 Min.   : 5.416  
 1st Qu.: 6.908  
 Median : 8.071  
 Mean   : 7.819  
 3rd Qu.: 9.486  
 Max.   :10.181  
                 

Code

summary(dt_prey)

   Study_ID          Cohort_ID         Shared_control_ID       Year     
 Length:146         Length:146         Length:146         Min.   :2003  
 Class :character   Class :character   Class :character   1st Qu.:2008  
 Mode  :character   Mode  :character   Mode  :character   Median :2012  
                                                          Mean   :2011  
                                                          3rd Qu.:2013  
                                                          Max.   :2022  
                                                                        
   Country          Data_source        Data_location      Bird_species      
 Length:146         Length:146         Length:146         Length:146        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
 Bird_common_name   Prey_species       Prey_common_name     Response        
 Length:146         Length:146         Length:146         Length:146        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
   Measure           Direction         Detail_result      Note_result       
 Length:146         Length:146         Length:146         Length:146        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
 Treatment_stimulus Number_pattern   Diameter_pattern  Area_pattern   
 Length:146         Min.   : 1.000   Min.   : 1.700   Min.   :  2.27  
 Class :character   1st Qu.: 2.000   1st Qu.: 4.500   1st Qu.: 15.90  
 Mode  :character   Median : 2.000   Median : 6.000   Median : 28.27  
                    Mean   : 2.397   Mean   : 6.741   Mean   : 36.04  
                    3rd Qu.: 2.000   3rd Qu.: 8.380   3rd Qu.: 43.01  
                    Max.   :11.000   Max.   :15.000   Max.   :176.71  
                                                                      
 Area_background    Area_ratio      Shape_pattern      Control_stimulus  
 Min.   : 129.6   Min.   :0.00816   Length:146         Length:146        
 1st Qu.: 625.0   1st Qu.:0.04661   Class :character   Class :character  
 Median : 931.8   Median :0.08501   Mode  :character   Mode  :character  
 Mean   : 890.1   Mean   :0.08681                                        
 3rd Qu.:1110.6   3rd Qu.:0.12315                                        
 Max.   :2784.7   Max.   :0.22893                                        
                                                                         
 Note_stimulus       Type_prey          Shape_prey         Diet_bird        
 Length:146         Length:146         Length:146         Length:146        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
   Bird_sex           Bird_age               Tn                Cn         
 Length:146         Length:146         Min.   :  7.708   Min.   :  7.708  
 Class :character   Class :character   1st Qu.: 20.000   1st Qu.: 20.000  
 Mode  :character   Mode  :character   Median : 80.000   Median : 80.000  
                                       Mean   : 65.610   Mean   : 65.616  
                                       3rd Qu.: 96.000   3rd Qu.: 96.000  
                                       Max.   :160.000   Max.   :160.000  
                                                                          
     T_mean           C_mean          MT_mean           MC_mean       
 Min.   :0.1148   Min.   :0.1552   Min.   :-6.1300   Min.   :-4.3900  
 1st Qu.:0.4656   1st Qu.:0.3357   1st Qu.:-1.2048   1st Qu.:-0.9296  
 Median :0.7374   Median :0.4496   Median :-0.5582   Median :-0.3664  
 Mean   :1.0794   Mean   :0.8256   Mean   :-0.6946   Mean   :-0.4626  
 3rd Qu.:1.2864   3rd Qu.:1.1137   3rd Qu.:-0.1148   3rd Qu.:-0.1552  
 Max.   :6.1300   Max.   :4.3900   Max.   : 0.8822   Max.   : 0.9584  
 NA's   :127      NA's   :127      NA's   :121       NA's   :121      
      T_sd             C_sd         T_proportion      C_proportion   
 Min.   :0.1403   Min.   :0.0864   Min.   :0.01598   Min.   :0.0100  
 1st Qu.:0.6219   1st Qu.:0.1422   1st Qu.:0.18705   1st Qu.:0.1100  
 Median :0.9192   Median :0.2532   Median :0.30000   Median :0.1621  
 Mean   :1.2667   Mean   :0.4709   Mean   :0.34403   Mean   :0.2490  
 3rd Qu.:1.8102   3rd Qu.:0.4555   3rd Qu.:0.45000   3rd Qu.:0.3608  
 Max.   :5.4821   Max.   :3.6547   Max.   :0.90000   Max.   :0.8333  
 NA's   :124      NA's   :124      NA's   :19        NA's   :19      
 Study_design         Dataset          Mean_median        T_mean_median   
 Length:146         Length:146         Length:146         Min.   :0.1085  
 Class :character   Class :character   Class :character   1st Qu.:0.4496  
 Mode  :character   Mode  :character   Mode  :character   Median :0.6934  
                                                          Mean   :1.0349  
                                                          3rd Qu.:1.2320  
                                                          Max.   :6.1300  
                                                          NA's   :127     
 C_mean_median    Type_of_variance_statistic      T_se            C_se       
 Min.   :0.1424   Length:146                 Min.   :0.040   Min.   :0.0200  
 1st Qu.:0.3137   Class :character           1st Qu.:0.055   1st Qu.:0.0350  
 Median :0.4112   Mode  :character           Median :0.065   Median :0.0450  
 Mean   :0.7883                              Mean   :0.185   Mean   :0.1225  
 3rd Qu.:1.0391                              3rd Qu.:0.195   3rd Qu.:0.1325  
 Max.   :4.3900                              Max.   :0.570   Max.   :0.3800  
 NA's   :127                                 NA's   :142     NA's   :142     
 T_variance_q1     T_variance_q3    C_variance_q1    C_variance_q3   
 Min.   :0.07075   Min.   :0.1651   Min.   :0.1136   Min.   :0.2096  
 1st Qu.:0.28640   1st Qu.:0.6584   1st Qu.:0.2594   1st Qu.:0.4316  
 Median :0.47294   Median :0.9294   Median :0.3159   Median :0.5600  
 Mean   :0.52460   Mean   :1.1200   Mean   :0.4762   Mean   :0.8185  
 3rd Qu.:0.72800   3rd Qu.:1.6352   3rd Qu.:0.7712   3rd Qu.:1.3267  
 Max.   :1.22170   Max.   :2.3216   Max.   :0.9481   Max.   :1.6880  
 NA's   :128       NA's   :128      NA's   :128      NA's   :128     
     Note                lnRR              lnRR_var            Obs_ID      
 Length:146         Min.   :-0.815486   Min.   :0.002201   Min.   :  1.00  
 Class :character   1st Qu.:-0.001665   1st Qu.:0.011064   1st Qu.: 37.25  
 Mode  :character   Median : 0.228398   Median :0.021493   Median : 73.50  
                    Mean   : 0.262303   Mean   :0.060098   Mean   : 73.50  
                    3rd Qu.: 0.479181   3rd Qu.:0.051982   3rd Qu.:109.75  
                    Max.   : 2.247699   Max.   :0.434523   Max.   :146.00  
                                                                           
  Log_diameter       Log_area      Log_background 
 Min.   :0.5306   Min.   :0.8197   Min.   :4.864  
 1st Qu.:1.5041   1st Qu.:2.7666   1st Qu.:6.438  
 Median :1.7918   Median :3.3420   Median :6.837  
 Mean   :1.8182   Mean   :3.3132   Mean   :6.667  
 3rd Qu.:2.1258   3rd Qu.:3.7614   3rd Qu.:7.013  
 Max.   :2.7081   Max.   :5.1745   Max.   :7.932  
                                                  

Code

summary(dt_all)

   Study_ID          Cohort_ID         Shared_control_ID       Year     
 Length:263         Length:263         Length:263         Min.   :1957  
 Class :character   Class :character   Class :character   1st Qu.:2009  
 Mode  :character   Mode  :character   Mode  :character   Median :2010  
                                                          Mean   :2009  
                                                          3rd Qu.:2013  
                                                          Max.   :2022  
                                                                        
   Country          Data_source        Data_location      Bird_species      
 Length:263         Length:263         Length:263         Length:263        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
 Bird_common_name   Prey_species       Prey_common_name     Response        
 Length:263         Length:263         Length:263         Length:263        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
   Measure           Direction         Detail_result      Note_result       
 Length:263         Length:263         Length:263         Length:263        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
 Treatment_stimulus Number_pattern   Diameter_pattern  Area_pattern   
 Length:263         Min.   : 1.000   Min.   : 1.70    Min.   :  2.27  
 Class :character   1st Qu.: 2.000   1st Qu.: 4.50    1st Qu.: 15.90  
 Mode  :character   Median : 2.000   Median : 7.00    Median : 38.48  
                    Mean   : 2.194   Mean   :10.55    Mean   :133.28  
                    3rd Qu.: 2.000   3rd Qu.:11.54    3rd Qu.: 81.50  
                    Max.   :11.000   Max.   :31.28    Max.   :490.87  
                                                                      
 Area_background     Area_ratio      Shape_pattern      Control_stimulus  
 Min.   :  129.6   Min.   :0.00816   Length:263         Length:263        
 1st Qu.:  697.5   1st Qu.:0.03872   Class :character   Class :character  
 Median : 1040.0   Median :0.07434   Mode  :character   Mode  :character  
 Mean   : 3097.0   Mean   :0.08947                                        
 3rd Qu.: 1444.0   3rd Qu.:0.12315                                        
 Max.   :26400.0   Max.   :0.33071                                        
                                                                          
 Note_stimulus       Type_prey          Shape_prey         Diet_bird        
 Length:263         Length:263         Length:263         Length:263        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
   Bird_sex           Bird_age               Tn               Cn        
 Length:263         Length:263         Min.   :  4.00   Min.   :  4.00  
 Class :character   Class :character   1st Qu.:  8.00   1st Qu.:  8.00  
 Mode  :character   Mode  :character   Median : 20.00   Median : 20.00  
                                       Mean   : 42.97   Mean   : 42.77  
                                       3rd Qu.: 82.00   3rd Qu.: 82.00  
                                       Max.   :160.00   Max.   :160.00  
                                                                        
     T_mean              C_mean             MT_mean         
 Min.   :   0.0004   Min.   :   0.0400   Min.   :-556.1000  
 1st Qu.:   0.6021   1st Qu.:   0.3675   1st Qu.:  -1.3776  
 Median :   7.2595   Median :   3.9950   Median :  -0.1148  
 Mean   : 100.4589   Mean   :  92.8595   Mean   : -26.1327  
 3rd Qu.:  93.1250   3rd Qu.:  56.6866   3rd Qu.:   0.8822  
 Max.   :1180.4878   Max.   :2115.0000   Max.   :   0.9996  
 NA's   :143         NA's   :143         NA's   :202        
    MC_mean                T_sd               C_sd           T_proportion   
 Min.   :-2115.0000   Min.   :  0.0100   Min.   :   0.050   Min.   :0.0100  
 1st Qu.:   -1.1749   1st Qu.:  0.5278   1st Qu.:   0.229   1st Qu.:0.1667  
 Median :   -0.1552   Median :  4.2319   Median :   2.711   Median :0.2984  
 Mean   :  -72.1179   Mean   : 67.7837   Mean   :  64.351   Mean   :0.3353  
 3rd Qu.:    0.8949   3rd Qu.: 59.4457   3rd Qu.:  52.998   3rd Qu.:0.4713  
 Max.   :    0.9584   Max.   :731.8482   Max.   :1052.175   Max.   :0.9000  
 NA's   :202          NA's   :140        NA's   :140        NA's   :94      
  C_proportion    Study_design         Dataset          Mean_median       
 Min.   :0.0100   Length:263         Length:263         Length:263        
 1st Qu.:0.1045   Class :character   Class :character   Class :character  
 Median :0.1770   Mode  :character   Mode  :character   Mode  :character  
 Mean   :0.2613                                                           
 3rd Qu.:0.4000                                                           
 Max.   :0.8333                                                           
 NA's   :94                                                               
 T_mean_median       C_mean_median       Type_of_variance_statistic
 Min.   :   0.0004   Min.   :   0.0000   Length:263                
 1st Qu.:   0.6584   1st Qu.:   0.4022   Class :character          
 Median :   7.2766   Median :   4.7365   Mode  :character          
 Mean   :  89.8277   Mean   :  83.4133                             
 3rd Qu.:  58.0122   3rd Qu.:  46.5061                             
 Max.   :1180.4878   Max.   :2115.0000                             
 NA's   :127         NA's   :127                                   
      T_se                C_se          T_variance_q1      T_variance_q3     
 Min.   :  0.00000   Min.   :  0.0000   Min.   : 0.07075   Min.   :  0.1651  
 1st Qu.:  0.08208   1st Qu.:  0.0877   1st Qu.: 0.28640   1st Qu.:  0.6656  
 Median :  3.52000   Median :  1.1000   Median : 0.50720   Median :  0.9296  
 Mean   : 24.82370   Mean   : 25.6820   Mean   : 4.57594   Mean   : 14.0085  
 3rd Qu.: 24.55107   3rd Qu.: 19.8747   3rd Qu.: 0.78560   3rd Qu.:  1.7799  
 Max.   :182.92683   Max.   :372.0000   Max.   :77.50000   Max.   :246.0000  
 NA's   :162         NA's   :162        NA's   :244        NA's   :244       
 C_variance_q1     C_variance_q3         Note                lnRR         
 Min.   : 0.1136   Min.   : 0.2096   Length:263         Min.   :-2.15769  
 1st Qu.: 0.2594   1st Qu.: 0.4340   Class :character   1st Qu.:-0.04159  
 Median : 0.3440   Median : 0.5936   Mode  :character   Median : 0.15315  
 Mean   : 1.9774   Mean   : 2.6702                      Mean   : 0.26794  
 3rd Qu.: 0.8048   3rd Qu.: 1.4151                      3rd Qu.: 0.49057  
 Max.   :29.0000   Max.   :36.0000                      Max.   : 4.16777  
 NA's   :244       NA's   :244                                            
    lnRR_var             Obs_ID       Log_diameter       Log_area     
 Min.   :0.0008295   Min.   :  1.0   Min.   :0.5306   Min.   :0.8197  
 1st Qu.:0.0108688   1st Qu.: 66.5   1st Qu.:1.5041   1st Qu.:2.7666  
 Median :0.0238924   Median :132.0   Median :1.9459   Median :3.6503  
 Mean   :0.0889528   Mean   :132.0   Mean   :2.0850   Mean   :3.8740  
 3rd Qu.:0.0743148   3rd Qu.:197.5   3rd Qu.:2.4456   3rd Qu.:4.4006  
 Max.   :1.3885984   Max.   :263.0   Max.   :3.4430   Max.   :6.1962  
                                                                      
 Log_background  
 Min.   : 4.864  
 1st Qu.: 6.542  
 Median : 6.947  
 Mean   : 7.179  
 3rd Qu.: 7.275  
 Max.   :10.181  
                 

Final dataset for meta-analysis

I cannot attach the caption to datatable() format table. I need to use kable() format table?

Code

datatable(dt_all, caption = "Dataset for meta-analysis", 
          options = list(pageLength = 10, scrollX = TRUE))

Table 2 - Dataset for meta-analysis

Check phylogenetic relatedness

If phylogeny do not explain heterogeniety much, we will not need to consider it and the two datasets can be merged.

We used meta-analytical models to calculate total I2 (a measure of heterogeneity not caused by sampling error; Higgins. et al., 2003) and the partial I2 explained by each random factor.

Based on this, we need not to consider phylogenetic relatedness from our meta-regressions as this random factor did not explain much of the heterogeneity between effect sizes (partial I2 < 0.001%), then we can combine two datasets (predator and prey) for meta-analysis.

Code

phy_model <- function(cor_tree = vcv_tree){
            model <- rma.mv(yi = lnRR,
                            V = VCV_pred,
                            random = list(~1 | Study_ID,
                                          ~1 | Shared_control_ID,
                                          ~1 | Cohort_ID,
                                          ~1 | Obs_ID,
                                          ~1 | Bird_species,
                                          ~1 | Phylogeny),
                            R = list(Phylogeny = cor_tree),
                            test = "t",
                            method = "REML",
                            sparse = TRUE,
                            data = dat_pred)
  model
}

tree_50 <- trees[1:50]
vcv_tree_50 <- map(tree_50, ~vcv(.x, corr = TRUE))

# Run multiple meta-analyses with 50 different trees and obtain the combined result

ma_50 <- mclapply(vcv_tree_50, phy_model, mc.cores = 8) 

Code

ma_50 <- readRDS(here("Rdata", "ma_50.RDS"))

# combining the results
est_50 <- map_dbl(ma_50, ~ .x$b[[1]])
se_50 <-  map_dbl(ma_50, ~ .x$se)
df_50 <- c(rbind(est_50, se_50))

# creating an array required by pool.mi
my_array <- array(df_50, dim = c(1, 2, 50))

com_res <- round(pool.mi(my_array), 4)

Code

knitr::kable(com_res, caption = "Combined result of 50 phylogenetic trees")

Combined result of 50 phylogenetic trees

pooled.mean	pooled.total.se	se.within	se.between	relative.inc.var	frac.miss.info	CI.1	CI.2	p.value
0.1394	0.1192	0.1192	0	0	0	-0.0943	0.3731	0.2424

Code

sigma2_mod <- do.call(rbind, lapply(ma_50, function(x) x$sigma2))
sigma2_mod <- data.frame(sigma2_mod)

colnames(sigma2_mod) <- c("sigma^2.1_Study_ID", "sigma^2.2_SharedControl_ID", 
                          "sigma^2.3_Cohort_ID", "sigma^2.4_Obs_ID", 
                          "sigma^2.5_BirdSpecies", "sigma^2.6_Phylogeny")

sigma2_mod <- round(colMeans(sigma2_mod), 4)

knitr::kable(sigma2_mod, caption = "Average variance components")

Average variance components

	x
sigma^2.1_Study_ID	0.0000
sigma^2.2_SharedControl_ID	0.0923
sigma^2.3_Cohort_ID	0.1009
sigma^2.4_Obs_ID	0.5323
sigma^2.5_BirdSpecies	0.0000
sigma^2.6_Phylogeny	0.0000

We got 1000 phylogenetic trees from BirdTree.org

Only 50 trees are used as in Nakagawa & Villemereuil (2019)

2 Meta-analysis - predator dataset

We used meta-analytical models to calculate total I2 (a measure of heterogeneity not caused by sampling error; Higgins. et al., 2003) and the partial I2 explained by each random factor.

Which random effects remove in our models?

We can remove Study_ID

Code

ma_pred <- rma.mv(yi = lnRR,
                  V = VCV_pred, 
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_pred)

summary(ma_pred)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-150.0443   300.0886   310.0886   323.8565   310.6340   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0000  0.0001     18     no           Study_ID 
sigma^2.2  0.0786  0.2803     30     no  Shared_control_ID 
sigma^2.3  0.1271  0.3565     34     no          Cohort_ID 
sigma^2.4  0.4563  0.6755    117     no             Obs_ID 

Test for Heterogeneity:
Q(df = 116) = 1280.1795, p-val < .0001

Model Results:

estimate      se    tval   df    pval    ci.lb   ci.ub    
  0.1657  0.1164  1.4234  116  0.1573  -0.0649  0.3963    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

pred_i2 <- round(i2_ml(ma_pred), 4)

pred_r2 <- round(r2_ml(ma_pred), 4)

I2_Total	I2_Study_ID	I2_Shared_control_ID	I2_Cohort_ID	I2_Obs_ID
98.8884	0	11.7364	18.9879	68.1641

R2_marginal	R2_conditional
0	0.3107

Table 3 - Heterogeneity of effect size

Code

orchard_plot(ma_pred,
              group = "Study_ID",
              xlab = "log response ratio (lnRR)", 
              angle = 45) + 
              scale_x_discrete(labels = c("Overall effect")) + 
              scale_colour_okabe_ito()+
              scale_fill_okabe_ito()

Figure 2— Estimates of overpred effect size and 95% confidence intervals

2.1 Meta-regressions: uni-moderator - predetor dataset

Pattern type

Eyespot vs. conspicuous patterns - Is there significant difference of effect size between two patterns?

Code

#normal model
pred_mr_eyespot <- rma.mv(yi = lnRR,
                    V = VCV_pred, 
                    mods = ~ Treatment_stimulus,
                    random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                    test = "t",
                    method = "REML", 
                    sparse = TRUE,
                    data = dt_pred)

summary(pred_mr_eyespot)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.9140   297.8280   307.8280   321.5527   308.3785   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0833  0.2885     30     no  Shared_control_ID 
sigma^2.2  0.1229  0.3506     34     no          Cohort_ID 
sigma^2.3  0.4600  0.6782    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1278.4557, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.6546, p-val = 0.4201

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.0809  0.3268  -0.2475  115  0.8049  -0.7283 
Treatment_stimuluseyespot    0.2812  0.3476   0.8091  115  0.4201  -0.4073 
                            ci.ub    
intrcpt                    0.5665    
Treatment_stimuluseyespot  0.9697    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

pred_r2_1 <- round(r2_ml(pred_mr_eyespot), 4)

R2_marginal	R2_conditional
0.0067	0.3141

Table XX. Model goodness-of-fit

Code

orchard_plot(pred_mr_eyespot,
            mod = "Treatment_stimulus",
            group = "Study_ID",
            xlab = "log response ratio (lnRR)",
            angle = 45) + 
            scale_colour_okabe_ito(order = 2:3)+
            scale_fill_okabe_ito(order = 2:3)

Figure 3— Eyespot vs. conspicuous patterns

Code

# intercept-removed model
pred_mr_eyespot1 <- rma.mv(yi = lnRR,
                      V = VCV_pred, 
                      mods = ~ Treatment_stimulus -1,
                      random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                                test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = dt_pred)

summary(pred_mr_eyespot1)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.9140   297.8280   307.8280   321.5527   308.3785   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0833  0.2885     30     no  Shared_control_ID 
sigma^2.2  0.1229  0.3506     34     no          Cohort_ID 
sigma^2.3  0.4600  0.6782    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1278.4557, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 115) = 1.3386, p-val = 0.2663

Model Results:

                               estimate      se     tval   df    pval    ci.lb 
Treatment_stimulusconspicuous   -0.0809  0.3268  -0.2475  115  0.8049  -0.7283 
Treatment_stimuluseyespot        0.2003  0.1242   1.6127  115  0.1096  -0.0457 
                                ci.ub    
Treatment_stimulusconspicuous  0.5665    
Treatment_stimuluseyespot      0.4463    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pattern diameter (mm)

Code

pred_mr_diameter <- rma.mv(yi = lnRR,
                      V = VCV_pred, 
                      mods = ~ Log_diameter,
                      random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = dt_pred)

summary(pred_mr_diameter)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.6526   297.3051   307.3051   321.0298   307.8556   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0829  0.2879     30     no  Shared_control_ID 
sigma^2.2  0.1316  0.3628     34     no          Cohort_ID 
sigma^2.3  0.4566  0.6758    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1269.1592, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.3781, p-val = 0.5398

Model Results:

              estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt        -0.0240  0.3296  -0.0729  115  0.9420  -0.6768  0.6288    
Log_diameter    0.0880  0.1432   0.6149  115  0.5398  -0.1955  0.3716    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

pred_r2_2 <- round(r2_ml(pred_mr_diameter), 4)

R2_marginal	R2_conditional
0.0085	0.3254

Table XX. Model goodness-of-fit

Code

bubble_plot(pred_mr_diameter,
            mod = "Log_diameter",
            group = "Study_ID",
            xlab = "Log-transformed diameter")

Pattern area (mm²)

Code

pred_mr_area <- rma.mv(yi = lnRR,
                  V = VCV_pred, 
                  mods = ~ Log_area,
                  random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_pred)

summary(pred_mr_area)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.5635   297.1269   307.1269   320.8516   307.6774   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0816  0.2857     30     no  Shared_control_ID 
sigma^2.2  0.1298  0.3603     34     no          Cohort_ID 
sigma^2.3  0.4571  0.6761    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1269.5060, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.5212, p-val = 0.4718

Model Results:

          estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt    -0.0476  0.3175  -0.1498  115  0.8812  -0.6766  0.5814    
Log_area    0.0529  0.0732   0.7219  115  0.4718  -0.0922  0.1979    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

pred_r2_3 <- round(r2_ml(pred_mr_area), 4)

R2_marginal	R2_conditional
0.0122	0.3245

Table XX. Model goodness-of-fit

Code

bubble_plot(pred_mr_area,
            mod = "Log_area",
            group = "Study_ID",
            xlab = "Log-transformed area")

Pattern number

Code

pred_mr_num <- rma.mv(yi = lnRR,
                  V = VCV_pred,
                  mods = ~ Number_pattern,
                  random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 |Obs_ID),                      
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_pred)

summary(pred_mr_num)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-149.1039   298.2077   308.2077   321.9324   308.7582   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0838  0.2894     30     no  Shared_control_ID 
sigma^2.2  0.1289  0.3590     34     no          Cohort_ID 
sigma^2.3  0.4602  0.6784    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1277.6943, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.0197, p-val = 0.8885

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt           0.2002  0.2728   0.7338  115  0.4645  -0.3402  0.7405    
Number_pattern   -0.0161  0.1145  -0.1405  115  0.8885  -0.2430  0.2108    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

pred_r2_4 <- round(r2_ml(pred_mr_num), 4)

R2_marginal	R2_conditional
3e-04	0.3162

Table XX. Model goodness-of-fit

Code

bubble_plot(pred_mr_num,
            mod = "Number_pattern",
            group = "Study_ID",
            xlab = "Number of patterns") 

Prey material type

Our dataset includes two material types of prey: live and artificial. Is there significant difference of effect size between two stimuli?

Code

# normal model
pred_mr_prey_type <- rma.mv(yi = lnRR,
                  V = VCV_pred, 
                  mods = ~ Type_prey,
                  random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID), 
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_pred)

summary(pred_mr_prey_type)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.7653   297.5306   307.5306   321.2553   308.0810   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0793  0.2816     30     no  Shared_control_ID 
sigma^2.2  0.1530  0.3912     34     no          Cohort_ID 
sigma^2.3  0.4487  0.6699    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1232.0173, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.4715, p-val = 0.4937

Model Results:

               estimate      se    tval   df    pval    ci.lb   ci.ub    
intrcpt          0.0994  0.1519  0.6543  115  0.5142  -0.2015  0.4002    
Type_preyreal    0.1699  0.2474  0.6866  115  0.4937  -0.3202  0.6600    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

pred_r2_5 <- round(r2_ml(pred_mr_prey_type), 4)

R2_marginal	R2_conditional
0.0079	0.3463

Table XX. Model goodness-of-fit

Code

orchard_plot(pred_mr_prey_type,
              mod = "Type_prey",
              group = "Study_ID",
              xlab = "Type of prey",
              angle = 45) +
            scale_colour_okabe_ito(order = 4:5)+
            scale_fill_okabe_ito(order = 4:5)

Code

# intercept-removed model
pred_mr_prey_type1 <- rma.mv(yi = lnRR,
                            V = VCV_pred, 
                            mods = ~ Type_prey -1,
                            random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID), 
                            test = "t",
                            method = "REML", 
                            sparse = TRUE,
                            data = dt_pred)

summary(pred_mr_prey_type1)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.7653   297.5306   307.5306   321.2553   308.0810   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0793  0.2816     30     no  Shared_control_ID 
sigma^2.2  0.1530  0.3912     34     no          Cohort_ID 
sigma^2.3  0.4487  0.6699    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1232.0173, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 115) = 1.1641, p-val = 0.3159

Model Results:

                     estimate      se    tval   df    pval    ci.lb   ci.ub    
Type_preyartificial    0.0994  0.1519  0.6543  115  0.5142  -0.2015  0.4002    
Type_preyreal          0.2693  0.1953  1.3785  115  0.1707  -0.1177  0.6562    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prey shape type

Our dataset includes 4 types of prey type: real butterfly, abstract butterfly, abstract caterpillar, and abstract prey. Is there significant difference of effect size between two stimuli?

Code

pred_mr_prey_shape <- rma.mv(yi = lnRR,
                        V = VCV_pred, 
                        mods = ~ Shape_prey -1,
                        random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),,
                                      test = "t",
                                      method = "REML", 
                                      sparse = TRUE,
                                      data = dt_pred)

summary(pred_mr_prey_shape)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-146.3709   292.7417   306.7417   325.8334   307.8084   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0998  0.3159     30     no  Shared_control_ID 
sigma^2.2  0.1719  0.4146     34     no          Cohort_ID 
sigma^2.3  0.4475  0.6689    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 113) = 1227.2585, p-val < .0001

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 113) = 0.5653, p-val = 0.6883

Model Results:

                                estimate      se    tval   df    pval    ci.lb 
Shape_preyabstract_butterfly      0.1090  0.3403  0.3203  113  0.7493  -0.5652 
Shape_preyabstract_caterpillar    0.0181  0.2757  0.0656  113  0.9478  -0.5280 
Shape_preyabstract_prey           0.1425  0.2407  0.5920  113  0.5551  -0.3344 
Shape_preybutterfly               0.2724  0.2028  1.3431  113  0.1819  -0.1294 
                                 ci.ub    
Shape_preyabstract_butterfly    0.7832    
Shape_preyabstract_caterpillar  0.5642    
Shape_preyabstract_prey         0.6194    
Shape_preybutterfly             0.6742    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

pred_r2_6 <- round(r2_ml(pred_mr_prey_shape), 4)

R2_marginal	R2_conditional
0.0088	0.3832

Table XX. Model goodness-of-fit

Code

orchard_plot(pred_mr_prey_shape,
              mod = "Shape_prey",
              group = "Study_ID",
              xlab = "Shape of prey",
              angle = 45) +
              scale_colour_okabe_ito(order = 6:9)+
              scale_fill_okabe_ito(order = 6:9)

Figure 4— Effect size and prey shape

2.2 Correlation visualisation and choose moderators for multi-moderator meta-regression in predator dataset

Before we run multi-moderator meta-regressions, we need to consider the correlation between moderators. Area, diameter and background seem to be correlated. We visualised the correlation between these variables.

Visualisation of correlation

Code

# correlation between continuous variables
pred_corr_cont <- round(cor(dt_pred[, c("Diameter_pattern", "Area_pattern", "Number_pattern", "Area_background")]),2)

pred_p1 <- ggcorrplot(pred_corr_cont, hc.order = TRUE, lab = TRUE, outline.col = "white", colors = c("#6D9EC1", "white", "#E46726"), title = "(a) Continuous variables")

pred_corr_cont_log <- round(cor(dt_pred[, c("Log_diameter", "Log_area", "Number_pattern", "Log_background")]),2)

pred_p2 <- ggcorrplot(pred_corr_cont_log, hc.order = TRUE, lab = TRUE, outline.col = "white", colors = c("#6D9EC1", "white", "#E46726"), title = "(b) Log-transormed continuous variables")

pred_p1 + pred_p2 +plot_layout(guides = 'collect')

Figure 5— Correlation between coninuous moderators

Code

summary(dt_pred$Treatment_stimulus)

   Length     Class      Mode 
      117 character character 

Code

pred_dat1 <- dt_pred %>%
  dplyr::select("Treatment_stimulus", "Type_prey", "Shape_prey")

pred_corr_cat <- GKtauDataframe(pred_dat1)
plot(pred_corr_cat)

Figure 6— Correlation between categorical variables

We should not include “Shape_prey” and “Type_prey” in the model at the same time. Because they are correlated.

Choose moderators

We used model R2 values to find better model and modelator VIF values to check multicollinearity between moderators. Higher R2 indicates better model and VIF > 2 indicates multicollinearity.

2.2.0.1 R2

Code

pred_r2_area <- rma.mv(yi = lnRR,
                  V = VCV_pred,
                  mods = ~ Log_area + Number_pattern,
                  random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML",
                  sparse = TRUE,
                  data = dt_pred)

pred_r2_diameter <- rma.mv(yi = lnRR,
                      V = VCV_pred,
                      mods = ~ Log_diameter + Number_pattern,
                      random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_pred)

# check r2 values
r2_ml(pred_r2_area)

   R2_marginal R2_conditional 
    0.01211581     0.33044840 

Code

r2_ml(pred_r2_diameter)

   R2_marginal R2_conditional 
   0.008825168    0.330370102 

It seems area is a better predictor than diameter .

2.3 Meta-regressions: multi-moderator

All potential moderators (full model)

Include all moderators in the model.

Code

pred_mr_all <- rma.mv(yi = lnRR,
                V = VCV_pred, 
                mods = ~ Treatment_stimulus + Log_area 
                + Number_pattern + Type_prey,
                random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_pred)

summary(pred_mr_all)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-144.9280   289.8560   305.8560   327.6040   307.2541   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0859  0.2931     30     no  Shared_control_ID 
sigma^2.2  0.1358  0.3685     34     no          Cohort_ID 
sigma^2.3  0.4645  0.6816    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 112) = 1225.0519, p-val < .0001

Test of Moderators (coefficients 2:5):
F(df1 = 4, df2 = 112) = 0.5654, p-val = 0.6883

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.3783  0.5230  -0.7234  112  0.4709  -1.4146 
Treatment_stimuluseyespot    0.2844  0.3530   0.8056  112  0.4222  -0.4151 
Log_area                     0.0815  0.0789   1.0329  112  0.3039  -0.0748 
Number_pattern              -0.0657  0.1246  -0.5271  112  0.5992  -0.3127 
Type_preyreal                0.2876  0.2785   1.0328  112  0.3039  -0.2642 
                            ci.ub    
intrcpt                    0.6579    
Treatment_stimuluseyespot  0.9839    
Log_area                   0.2378    
Number_pattern             0.1813    
Type_preyreal              0.8395    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(pred_mr_all)

   R2_marginal R2_conditional 
    0.03194282     0.34468854 

Moderators associated with conspicuousness

Only include moderators related to conspicuousness (pattern area and number of patterns) in the model.

Code

pred_mr_cons <- rma.mv(yi = lnRR,
                V = VCV_pred, 
                mods = ~ Log_area + Number_pattern,
                random = list(~1 | Shared_control_ID,
                              ~1 | Cohort_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_pred)

summary(pred_mr_cons)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-147.6172   295.2344   307.2344   323.6516   308.0195   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0865  0.2942     30     no  Shared_control_ID 
sigma^2.2  0.1326  0.3642     34     no          Cohort_ID 
sigma^2.3  0.4609  0.6789    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 114) = 1268.2068, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 114) = 0.2623, p-val = 0.7698

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt          -0.0122  0.4025  -0.0304  114  0.9758  -0.8096  0.7852    
Log_area          0.0529  0.0742   0.7126  114  0.4775  -0.0941  0.1998    
Number_pattern   -0.0164  0.1152  -0.1425  114  0.8869  -0.2447  0.2118    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(pred_mr_cons)

   R2_marginal R2_conditional 
    0.01211581     0.33044840 

Code

bubble_plot(pred_mr_cons,
            mod = "Log_area",
            group = "Study_ID",
            xlab = "Log-transformed area")

Code

pred_mr_cons_1 <- rma.mv(yi = lnRR,
                V = VCV_pred, 
                mods = ~ Log_area + Number_pattern + Treatment_stimulus,
                random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_pred)

summary(pred_mr_cons_1)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-146.4650   292.9300   306.9300   326.0217   307.9967   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0903  0.3005     30     no  Shared_control_ID 
sigma^2.2  0.1320  0.3633     34     no          Cohort_ID 
sigma^2.3  0.4635  0.6808    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 113) = 1264.4653, p-val < .0001

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 113) = 0.3975, p-val = 0.7550

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.2715  0.5129  -0.5293  113  0.5976  -1.2877 
Log_area                     0.0548  0.0746   0.7340  113  0.4644  -0.0931 
Number_pattern              -0.0177  0.1157  -0.1530  113  0.8787  -0.2470 
Treatment_stimuluseyespot    0.2902  0.3529   0.8222  113  0.4127  -0.4090 
                            ci.ub    
intrcpt                    0.7447    
Log_area                   0.2027    
Number_pattern             0.2116    
Treatment_stimuluseyespot  0.9894    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(pred_mr_cons_1)

   R2_marginal R2_conditional 
    0.02106248     0.33834441 

Moderators associated with prey in predator dataset

Code

pred_mr_pre <- rma.mv(yi = lnRR,
                V = VCV_pred, 
                mods = ~ Treatment_stimulus + Type_prey,
                random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_pred)

summary(pred_mr_pre)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-147.6573   295.3146   307.3146   323.7318   308.0997   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0844  0.2905     30     no  Shared_control_ID 
sigma^2.2  0.1473  0.3838     34     no          Cohort_ID 
sigma^2.3  0.4534  0.6733    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 114) = 1230.6103, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 114) = 0.5265, p-val = 0.5921

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.1347  0.3416  -0.3943  114  0.6941  -0.8115 
Treatment_stimuluseyespot    0.2714  0.3530   0.7689  114  0.4435  -0.4278 
Type_preyreal                0.1622  0.2481   0.6536  114  0.5147  -0.3294 
                            ci.ub    
intrcpt                    0.5421    
Treatment_stimuluseyespot  0.9706    
Type_preyreal              0.6538    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(pred_mr_pre)

   R2_marginal R2_conditional 
    0.01177801     0.34603083 

2.4 Publication bias

Funnel plot

Code

# funnel plot - standard error
funnel(ma_pred, yaxis = "sei", 
      xlab = "Standarised residuals",
      ylab = "Precision (inverse of SE)" )

# funnel plot - inverse of standard error
pred_f2 <- funnel(ma_pred, yaxis = "seinv", 
      xlab = "Standarized residuals",
      ylab = "Precision (inverse of SE)",  col = c(alpha("orange", 0.5)))

Figure 7— Effect size and its standard error

Figure 8— Effect size and its inverse standard error

Figure 9— Funnel plot of effect size and its standard error

Egger’s test

Code

pred_df_bias <- dt_pred %>% mutate(sqrt_inv_e_n = sqrt((Cn + Tn)/(Cn * Tn)))

pred_bias_model <- rma.mv(yi = lnRR,
                      V = lnRR_var, 
                      mods = ~ 1 + sqrt_inv_e_n,
                      random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = pred_df_bias)

summary(pred_bias_model)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-149.0844   298.1687   308.1687   321.8934   308.7192   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0804  0.2835     30     no  Shared_control_ID 
sigma^2.2  0.1285  0.3585     34     no          Cohort_ID 
sigma^2.3  0.4601  0.6783    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1249.3192, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.0592, p-val = 0.8082

Model Results:

              estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt         0.2525  0.3752   0.6729  115  0.5023  -0.4907  0.9957    
sqrt_inv_e_n   -0.2161  0.8884  -0.2433  115  0.8082  -1.9758  1.5435    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(pred_bias_model,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size")

Figure 10— Egger’s test of publication bias

Decline effect

Code

pred_year_model <- rma.mv(yi = lnRR,
                      V = VCV_pred, 
                      mods = ~  1 + Year,
                      random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = pred_df_bias)

summary(pred_year_model)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.7939   297.5878   307.5878   321.3125   308.1383   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0843  0.2904     30     no  Shared_control_ID 
sigma^2.2  0.1359  0.3687     34     no          Cohort_ID 
sigma^2.3  0.4563  0.6755    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1254.3394, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.0430, p-val = 0.8361

Model Results:

         estimate       se     tval   df    pval     ci.lb    ci.ub    
intrcpt    4.3742  20.3031   0.2154  115  0.8298  -35.8423  44.5907    
Year      -0.0021   0.0101  -0.2073  115  0.8361   -0.0221   0.0179    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(pred_year_model,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication")

Figure 11— Decline effect

multi-moderator model

Code

pred_bias_all <- rma.mv(yi = lnRR,
                      V = VCV_pred, 
                      mods = ~ 1 + sqrt_inv_e_n + Year,
                      random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = pred_df_bias)

summary(pred_bias_all)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-147.7933   295.5867   307.5867   324.0038   308.3717   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0852  0.2919     30     no  Shared_control_ID 
sigma^2.2  0.1372  0.3704     34     no          Cohort_ID 
sigma^2.3  0.4606  0.6786    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 114) = 1225.6345, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 114) = 0.0569, p-val = 0.9447

Model Results:

              estimate       se     tval   df    pval     ci.lb    ci.ub    
intrcpt         5.6563  20.9515   0.2700  114  0.7877  -35.8483  47.1610    
sqrt_inv_e_n   -0.2454   0.9196  -0.2668  114  0.7901   -2.0670   1.5763    
Year           -0.0027   0.0104  -0.2584  114  0.7966   -0.0233   0.0179    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.5 Additional analyses

2.5.1 Background area (mm²)

Code

pred_mr_background <- rma.mv(yi = lnRR,
                        V = VCV_pred, 
                        mods = ~ Log_background,
                  random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                        test = "t",
                        method = "REML",
                        sparse = TRUE,
                        data = dt_pred)

summary(pred_mr_background)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.6517   297.3033   307.3033   321.0280   307.8538   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0960  0.3098     30     no  Shared_control_ID 
sigma^2.2  0.1420  0.3768     34     no          Cohort_ID 
sigma^2.3  0.4490  0.6701    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1270.4543, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.2344, p-val = 0.6292

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt           0.4975  0.6997   0.7110  115  0.4785  -0.8885  1.8834    
Log_background   -0.0446  0.0922  -0.4842  115  0.6292  -0.2272  0.1380    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(pred_mr_background,
            mod = "Log_background",
            group = "Study_ID",
            xlab = "Log-transformed backgroud")

Figure 12— Effect size and background area

2.6 Figure galary

Discrete moderators

Code

# overall effect
pred_p1 <- orchard_plot(ma_pred,
                    group = "Study_ID",
                    xlab = "log response ratio (lnRR)", 
                    angle = 45,
                    legend.pos = "bottom.right") + 
                    scale_x_discrete(labels = c("Overall effect")) + 
                    theme(legend.direction = "horizontal",
                    legend.text = element_text(size = 6)) +
                    
                    scale_colour_okabe_ito()+
                    scale_fill_okabe_ito()

#eyespot vs. conspicuous patterns
pred_p2 <- orchard_plot(pred_mr_eyespot,
                  mod = "Treatment_stimulus",
                  group = "Study_ID",
                  xlab = "log response ratio (lnRR)",
                  angle = 45,
                  legend.pos = "bottom.right") +
                  theme(legend.direction = "horizontal",
                  legend.text = element_text(size = 6)) +
                  scale_colour_okabe_ito(order = 2:3)+
                  scale_fill_okabe_ito(order = 2:3)
                  
# prey type - real vs. artificial
pred_p3 <- orchard_plot(pred_mr_prey_type,
                  mod = "Type_prey",
                  group = "Study_ID",
                  xlab = "log response ratio (lnRR)",
                  angle = 45,
                  legend.pos = "bottom.right") +
                  theme(legend.direction = "horizontal",
                  legend.text = element_text(size = 6))+ 
                  scale_colour_okabe_ito(order = 4:5) +
                  scale_fill_okabe_ito(order = 4:5)

# combine
pred_patch1 <- pred_p1 | (pred_p2 / pred_p3)
pred_patch1 + plot_annotation(tag_levels = "a")

# output figure as pdf
# ggsave("fig1.pdf", width = 10, height = 5, dpi = 450)

Figure 13— Discrete moderators

Continuous moderators

Code

# these figs were created by multi meta-regression model results
# log-transformed area
pred_p4 <- bubble_plot(pred_mr_area,
                  mod = "Log_area",
                  group = "Study_ID",
                  xlab = "Log-transformed area") +
                  scale_x_continuous(breaks = seq(0,7,1.5)) +
                  scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5))

# number of patterns
pred_p5 <- bubble_plot(pred_mr_num,
            mod = "Number_pattern",
            group = "Study_ID",
            xlab = "Number of patterns") +
                  scale_x_continuous(breaks = seq(0,11,1.5)) +
                  scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5))

# combine
pred_p4 / pred_p5 + plot_annotation(tag_levels = "a")

# output figure as pdf
# ggsave("fig2_multi.pdf", width = 10, height = 10, dpi = 450)

Figure 14— Continuous moderators

Publication bias

Code

# funnel plot

# pdf("fig3.pdf")

funnel(ma_pred, yaxis = "seinv", 
      xlab = "Standarised residuals",
      ylab = "Precision (inverse of SE)",
      xlim = c(-4.0, 4.5), 
      ylim = c(0.01, 80.0), 
      col = c(alpha("mediumvioletred", 0.5)),
      cex = 0.7)

# dev.off()

# Egger's test and decline effect
pred_p7 <- bubble_plot(pred_bias_model,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size") +
            scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5)) 

pred_p8 <- bubble_plot(pred_year_model,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication") +  scale_y_continuous(breaks = seq(-2.5, 4.0, 1.5)) 

# combine plots created by bubble_plot()
pred_pub <- pred_p7 / pred_p8
pred_pub +  plot_annotation(tag_levels = 'a') 

# output figure as pdf
# ggsave("fig4.pdf", width = 10, height = 10, dpi = 450)

Figure 15— Funnel plot

Figure 16— Egger’s test and Decline effect

Figure 17— Publication bias

3 Meta-analysis - prey dataset

We used meta-analytical models to calculate total I2 (a measure of heterogeneity not caused by sampling error; Higgins. et al., 2003) and the partial I2 explained by each random factor.

Code

ma_prey <- rma.mv(yi = lnRR,
                  V = VCV_prey, 
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_prey)

summary(ma_prey)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-53.1702  106.3404  116.3404  131.2241  116.7721   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0842  0.2901     15     no           Study_ID 
sigma^2.2  0.0209  0.1446     58     no  Shared_control_ID 
sigma^2.3  0.0054  0.0733    123     no          Cohort_ID 
sigma^2.4  0.0323  0.1798    146     no             Obs_ID 

Test for Heterogeneity:
Q(df = 145) = 897.2723, p-val < .0001

Model Results:

estimate      se    tval   df    pval   ci.lb   ci.ub     
  0.2335  0.0865  2.6990  145  0.0078  0.0625  0.4045  ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

prey_i2 <- round(i2_ml(ma_prey), 4)

prey_r2 <- round(r2_ml(ma_prey), 4)

I2_Total	I2_Study_ID	I2_Shared_control_ID	I2_Cohort_ID	I2_Obs_ID
90.8149	53.531	13.299	3.417	20.5678

R2_marginal	R2_conditional
0	0.7735

Table 3 - Heterogeneity of effect size

Code

orchard_plot(ma_prey,
              group = "Study_ID",
              xlab = "log response ratio (lnRR)", 
              angle = 45) + 
              scale_x_discrete(labels = c("Overall effect")) + 
              scale_colour_okabe_ito()+
              scale_fill_okabe_ito()

Figure 18— Estimates of overprey effect size and 95% confidence intervals

3.1 Meta-regressions: uni-moderator - preyetor dataset

Pattern type

Eyespot vs. conspicuous patterns - Is there significant difference of effect size between two patterns?

Code

#normal model
prey_mr_eyespot <- rma.mv(yi = lnRR,
                    V = VCV_prey, 
                    mods = ~ Treatment_stimulus,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                    test = "t",
                    method = "REML", 
                    sparse = TRUE,
                    data = dt_prey)

summary(prey_mr_eyespot)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-52.3565  104.7130  116.7130  134.5319  117.3261   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0916  0.3026     15     no           Study_ID 
sigma^2.2  0.0243  0.1558     58     no  Shared_control_ID 
sigma^2.3  0.0036  0.0601    123     no          Cohort_ID 
sigma^2.4  0.0316  0.1776    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 144) = 890.1002, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 144) = 1.4023, p-val = 0.2383

Model Results:

                           estimate      se    tval   df    pval    ci.lb 
intrcpt                      0.1793  0.1011  1.7736  144  0.0782  -0.0205 
Treatment_stimuluseyespot    0.0922  0.0779  1.1842  144  0.2383  -0.0617 
                            ci.ub    
intrcpt                    0.3790  . 
Treatment_stimuluseyespot  0.2461    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

prey_r2_1 <- round(r2_ml(prey_mr_eyespot), 4)

R2_marginal	R2_conditional
0.0139	0.7939

Table XX. Model goodness-of-fit

Code

orchard_plot(prey_mr_eyespot,
            mod = "Treatment_stimulus",
            group = "Study_ID",
            xlab = "log response ratio (lnRR)",
            angle = 45) + 
            scale_colour_okabe_ito(order = 2:3)+
            scale_fill_okabe_ito(order = 2:3)

Figure 19— Eyespot vs. conspicuous patterns

Code

# intercept-removed model
prey_mr_eyespot1 <- rma.mv(yi = lnRR,
                      V = VCV_prey, 
                      mods = ~ Treatment_stimulus -1,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = dt_prey)

summary(prey_mr_eyespot1)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-52.3565  104.7130  116.7130  134.5319  117.3261   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0916  0.3026     15     no           Study_ID 
sigma^2.2  0.0243  0.1558     58     no  Shared_control_ID 
sigma^2.3  0.0036  0.0601    123     no          Cohort_ID 
sigma^2.4  0.0316  0.1776    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 144) = 890.1002, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 144) = 4.1054, p-val = 0.0185

Model Results:

                               estimate      se    tval   df    pval    ci.lb 
Treatment_stimulusconspicuous    0.1793  0.1011  1.7736  144  0.0782  -0.0205 
Treatment_stimuluseyespot        0.2715  0.0951  2.8540  144  0.0050   0.0835 
                                ci.ub     
Treatment_stimulusconspicuous  0.3790   . 
Treatment_stimuluseyespot      0.4595  ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pattern diameter (mm)

Code

prey_mr_diameter <- rma.mv(yi = lnRR,
                      V = VCV_prey, 
                      mods = ~ Log_diameter,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = dt_prey)

summary(prey_mr_diameter)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-39.1909   78.3818   90.3818  108.2007   90.9949   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0280  0.1674     15     no           Study_ID 
sigma^2.2  0.0157  0.1252     58     no  Shared_control_ID 
sigma^2.3  0.0000  0.0000    123     no          Cohort_ID 
sigma^2.4  0.0318  0.1783    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 144) = 627.4798, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 144) = 33.8104, p-val < .0001

Model Results:

              estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt        -0.5391  0.1413  -3.8153  144  0.0002  -0.8185  -0.2598  *** 
Log_diameter    0.4314  0.0742   5.8147  144  <.0001   0.2848   0.5781  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

prey_r2_2 <- round(r2_ml(prey_mr_diameter), 4)

R2_marginal	R2_conditional
0.3334	0.7192

Table XX. Model goodness-of-fit

Code

bubble_plot(prey_mr_diameter,
            mod = "Log_diameter",
            group = "Study_ID",
            xlab = "Log-transformed diameter")

Pattern area (mm²)

Code

prey_mr_area <- rma.mv(yi = lnRR,
                  V = VCV_prey, 
                  mods = ~ Log_area,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_prey)

summary(prey_mr_area)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-35.4140   70.8280   82.8280  100.6469   83.4411   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0331  0.1820     15     no           Study_ID 
sigma^2.2  0.0150  0.1225     58     no  Shared_control_ID 
sigma^2.3  0.0000  0.0000    123     no          Cohort_ID 
sigma^2.4  0.0283  0.1682    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 144) = 632.4105, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 144) = 41.4866, p-val < .0001

Model Results:

          estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt    -0.5926  0.1388  -4.2686  144  <.0001  -0.8671  -0.3182  *** 
Log_area    0.2557  0.0397   6.4410  144  <.0001   0.1772   0.3341  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

prey_r2_3 <- round(r2_ml(prey_mr_area), 4)

R2_marginal	R2_conditional
0.3727	0.7677

Table XX. Model goodness-of-fit

Code

bubble_plot(prey_mr_area,
            mod = "Log_area",
            group = "Study_ID",
            xlab = "Log-transformed area")

Pattern number

Code

prey_mr_num <- rma.mv(yi = lnRR,
                  V = VCV_prey,
                  mods = ~ Number_pattern,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_prey)

summary(prey_mr_num)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-42.6511   85.3022   97.3022  115.1211   97.9154   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0903  0.3004     15     no           Study_ID 
sigma^2.2  0.0110  0.1049     58     no  Shared_control_ID 
sigma^2.3  0.0000  0.0000    123     no          Cohort_ID 
sigma^2.4  0.0299  0.1730    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 144) = 826.5550, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 144) = 24.3477, p-val < .0001

Model Results:

                estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt           0.4795  0.0999   4.7982  144  <.0001   0.2820   0.6770  *** 
Number_pattern   -0.0795  0.0161  -4.9343  144  <.0001  -0.1113  -0.0476  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

prey_r2_4 <- round(r2_ml(prey_mr_num), 4)

R2_marginal	R2_conditional
0.1055	0.7959

Table XX. Model goodness-of-fit

Code

bubble_plot(prey_mr_num,
            mod = "Number_pattern",
            group = "Study_ID",
            xlab = "Number of patterns") 

Prey material type

Our dataset includes two material types of prey: live and artificial. Is there significant difference of effect size between two stimuli?

Code

# normal model
prey_mr_prey_type <- rma.mv(yi = lnRR,
                  V = VCV_prey, 
                  mods = ~ Type_prey,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_prey)

summary(prey_mr_prey_type)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-52.1781  104.3562  116.3562  134.1751  116.9694   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0908  0.3013     15     no           Study_ID 
sigma^2.2  0.0212  0.1455     58     no  Shared_control_ID 
sigma^2.3  0.0049  0.0699    123     no          Cohort_ID 
sigma^2.4  0.0326  0.1805    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 144) = 894.6573, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 144) = 0.3374, p-val = 0.5623

Model Results:

               estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt          0.2719  0.1099   2.4741  144  0.0145   0.0547  0.4892  * 
Type_preyreal   -0.1092  0.1881  -0.5808  144  0.5623  -0.4809  0.2625    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

prey_r2_5 <- round(r2_ml(prey_mr_prey_type), 4)

R2_marginal	R2_conditional
0.013	0.7848

Table XX. Model goodness-of-fit

Code

orchard_plot(prey_mr_prey_type,
              mod = "Type_prey",
              group = "Study_ID",
              xlab = "Type of prey",
              angle = 45) +
            scale_colour_okabe_ito(order = 4:5)+
            scale_fill_okabe_ito(order = 4:5)

Code

# intercept-removed model
prey_mr_prey_type1 <- rma.mv(yi = lnRR,
                            V = VCV_prey, 
                            mods = ~ Type_prey -1,
                            random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                            test = "t",
                            method = "REML", 
                            sparse = TRUE,
                            data = dt_prey)

summary(prey_mr_prey_type1)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-52.1781  104.3562  116.3562  134.1751  116.9694   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0908  0.3013     15     no           Study_ID 
sigma^2.2  0.0212  0.1455     58     no  Shared_control_ID 
sigma^2.3  0.0049  0.0699    123     no          Cohort_ID 
sigma^2.4  0.0326  0.1805    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 144) = 894.6573, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 144) = 3.6288, p-val = 0.0290

Model Results:

                     estimate      se    tval   df    pval    ci.lb   ci.ub    
Type_preyartificial    0.2719  0.1099  2.4741  144  0.0145   0.0547  0.4892  * 
Type_preyreal          0.1627  0.1526  1.0661  144  0.2882  -0.1389  0.4643    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prey shape type

Our dataset includes 4 types of prey type: real butterfly, abstract butterfly, abstract caterpillar, and abstract prey. Is there significant difference of effect size between two stimuli?

Code

prey_mr_prey_shape <- rma.mv(yi = lnRR,
                        V = VCV_prey, 
                        mods = ~ Shape_prey -1,
                        random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                        test = "t",
                        method = "REML", 
                        sparse = TRUE,
                        data = dt_prey)

summary(prey_mr_prey_shape)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-49.5600   99.1201  113.1201  133.8600  113.9497   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0701  0.2647     15     no           Study_ID 
sigma^2.2  0.0224  0.1498     58     no  Shared_control_ID 
sigma^2.3  0.0045  0.0671    123     no          Cohort_ID 
sigma^2.4  0.0324  0.1800    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 143) = 765.8269, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 143) = 3.9811, p-val = 0.0093

Model Results:

                                estimate      se    tval   df    pval    ci.lb 
Shape_preyabstract_butterfly      0.4043  0.1238  3.2662  143  0.0014   0.1596 
Shape_preyabstract_caterpillar    0.0252  0.1639  0.1538  143  0.8780  -0.2988 
Shape_preybutterfly               0.1568  0.1402  1.1187  143  0.2652  -0.1203 
                                 ci.ub     
Shape_preyabstract_butterfly    0.6490  ** 
Shape_preyabstract_caterpillar  0.3493     
Shape_preybutterfly             0.4339     

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

prey_r2_6 <- round(r2_ml(prey_mr_prey_shape), 4)

R2_marginal	R2_conditional
0.1853	0.796

Table XX. Model goodness-of-fit

Code

orchard_plot(prey_mr_prey_shape,
              mod = "Shape_prey",
              group = "Study_ID",
              xlab = "Shape of prey",
              angle = 45) +
              scale_colour_okabe_ito(order = 6:9)+
              scale_fill_okabe_ito(order = 6:9)

Figure 20— Effect size and prey shape

3.2 Correlation visualisation and choose moderators for multi-moderator meta-regression in preyator dataset

Before we run multi-moderator meta-regressions, we need to consider the correlation between moderators. Area, diameter and background seem to be correlated. We visualised the correlation between these variables.

Visualisation of correlation

Code

# correlation between continuous variables
prey_corr_cont <- round(cor(dt_prey[, c("Diameter_pattern", "Area_pattern", "Number_pattern", "Area_background")]),2)

prey_p1 <- ggcorrplot(prey_corr_cont, hc.order = TRUE, lab = TRUE, outline.col = "white", colors = c("#6D9EC1", "white", "#E46726"), title = "(a) Continuous variables")

prey_corr_cont_log <- round(cor(dt_prey[, c("Log_diameter", "Log_area", "Number_pattern", "Log_background")]),2)

prey_p2 <- ggcorrplot(prey_corr_cont_log, hc.order = TRUE, lab = TRUE, outline.col = "white", colors = c("#6D9EC1", "white", "#E46726"), title = "(b) Log-transormed continuous variables")

prey_p1 + prey_p2 +plot_layout(guides = 'collect')

Figure 21— Correlation between coninuous moderators

Code

prey_dat1 <- dt_prey %>%
  dplyr::select("Treatment_stimulus", "Type_prey", "Shape_prey")

prey_corr_cat <- GKtauDataframe(prey_dat1)
plot(prey_corr_cat)

Figure 22— Correlation between categorical variables

We should not include “Shape_prey” and “Type_prey” in the model at the same time. Because they are correlated.

Choose moderators

We used model R2 values to find better model and modelator VIF values to check multicollinearity between moderators. Higher R2 indicates better model and VIF > 2 indicates multicollinearity.

3.2.0.1 R2

Code

prey_r2_area <- rma.mv(yi = lnRR,
                  V = VCV_prey,
                  mods = ~ Log_area + Number_pattern,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML",
                  sparse = TRUE,
                  data = dt_prey)

prey_r2_diameter <- rma.mv(yi = lnRR,
                      V = VCV_prey,
                      mods = ~ Log_diameter + Number_pattern,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_prey)

# check r2 values
r2_ml(prey_r2_area)

   R2_marginal R2_conditional 
     0.3262005      0.7828015 

Code

r2_ml(prey_r2_diameter)

   R2_marginal R2_conditional 
     0.2741172      0.7602716 

It seems area is a better preyictor than diameter .

3.3 Meta-regressions: multi-moderator

All potential moderators (full model)

Include all moderators in the model.

Code

prey_mr_all <- rma.mv(yi = lnRR,
                V = VCV_prey, 
                mods = ~ Treatment_stimulus + Log_area 
                + Number_pattern + Type_prey,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_prey)

summary(prey_mr_all)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-32.5310   65.0620   83.0620  109.6009   84.4361   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0537  0.2318     15     no           Study_ID 
sigma^2.2  0.0132  0.1150     58     no  Shared_control_ID 
sigma^2.3  0.0000  0.0000    123     no          Cohort_ID 
sigma^2.4  0.0254  0.1592    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 141) = 603.8687, p-val < .0001

Test of Moderators (coefficients 2:5):
F(df1 = 4, df2 = 141) = 11.5068, p-val < .0001

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.4005  0.2099  -1.9080  141  0.0584  -0.8155 
Treatment_stimuluseyespot    0.0972  0.0706   1.3769  141  0.1707  -0.0424 
Log_area                     0.2082  0.0493   4.2273  141  <.0001   0.1108 
Number_pattern              -0.0384  0.0192  -2.0038  141  0.0470  -0.0763 
Type_preyreal                0.0851  0.1630   0.5219  141  0.6026  -0.2372 
                             ci.ub      
intrcpt                     0.0145    . 
Treatment_stimuluseyespot   0.2368      
Log_area                    0.3056  *** 
Number_pattern             -0.0005    * 
Type_preyreal               0.4073      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(prey_mr_all)

   R2_marginal R2_conditional 
     0.2973487      0.8069819 

Moderators associated with conspicuousness

Only include moderators related to conspicuousness (pattern area and number of patterns) in the model.

Code

prey_mr_cons <- rma.mv(yi = lnRR,
                V = VCV_prey, 
                mods = ~ Log_area + Number_pattern,
                random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_prey)

summary(prey_mr_cons)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-34.2630   68.5260   82.5260  103.2659   83.3556   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0454  0.2131     15     no           Study_ID 
sigma^2.2  0.0120  0.1095     58     no  Shared_control_ID 
sigma^2.3  0.0000  0.0000    123     no          Cohort_ID 
sigma^2.4  0.0273  0.1653    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 143) = 626.5275, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 143) = 21.9084, p-val < .0001

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub      
intrcpt          -0.3378  0.1999  -1.6893  143  0.0933  -0.7330  0.0575    . 
Log_area          0.2084  0.0484   4.3085  143  <.0001   0.1128  0.3040  *** 
Number_pattern   -0.0328  0.0185  -1.7769  143  0.0777  -0.0694  0.0037    . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(prey_mr_cons)

   R2_marginal R2_conditional 
     0.3262005      0.7828015 

Code

bubble_plot(prey_mr_cons,
            mod = "Log_area",
            group = "Study_ID",
            xlab = "Log-transformed area")

Code

prey_mr_cons_1 <- rma.mv(yi = lnRR,
                V = VCV_prey, 
                mods = ~ Log_area + Number_pattern + Treatment_stimulus,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_prey)

summary(prey_mr_cons_1)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-33.0836   66.1672   82.1672  105.8138   83.2499   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0464  0.2154     15     no           Study_ID 
sigma^2.2  0.0137  0.1169     58     no  Shared_control_ID 
sigma^2.3  0.0000  0.0000    123     no          Cohort_ID 
sigma^2.4  0.0253  0.1592    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 142) = 604.9523, p-val < .0001

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 142) = 15.4722, p-val < .0001

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.3948  0.2037  -1.9381  142  0.0546  -0.7975 
Log_area                     0.2089  0.0485   4.3035  142  <.0001   0.1129 
Number_pattern              -0.0349  0.0185  -1.8828  142  0.0618  -0.0716 
Treatment_stimuluseyespot    0.1065  0.0674   1.5798  142  0.1164  -0.0268 
                            ci.ub      
intrcpt                    0.0079    . 
Log_area                   0.3049  *** 
Number_pattern             0.0017    . 
Treatment_stimuluseyespot  0.2398      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(prey_mr_cons_1)

   R2_marginal R2_conditional 
     0.3171656      0.7974200 

Moderators associated with prey in preyator dataset

Code

prey_mr_pre <- rma.mv(yi = lnRR,
                V = VCV_prey, 
                mods = ~ Treatment_stimulus + Type_prey,
                random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_prey)

summary(prey_mr_pre)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-51.1934  102.3869  116.3869  137.1268  117.2165   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0940  0.3066     15     no           Study_ID 
sigma^2.2  0.0256  0.1599     58     no  Shared_control_ID 
sigma^2.3  0.0025  0.0501    123     no          Cohort_ID 
sigma^2.4  0.0318  0.1784    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 143) = 881.4144, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 143) = 1.1617, p-val = 0.3159

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                      0.2294  0.1163   1.9723  143  0.0505  -0.0005 
Treatment_stimuluseyespot    0.1127  0.0801   1.4065  143  0.1617  -0.0457 
Type_preyreal               -0.1816  0.1976  -0.9189  143  0.3597  -0.5721 
                            ci.ub    
intrcpt                    0.4593  . 
Treatment_stimuluseyespot  0.2711    
Type_preyreal              0.2090    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(prey_mr_pre)

   R2_marginal R2_conditional 
     0.0287115      0.7991167 

3.4 Publication bias

Funnel plot

Code

# funnel plot - standard error
funnel(ma_prey, yaxis = "sei", 
      xlab = "Standarised residuals",
      ylab = "Precision (inverse of SE)" )

# funnel plot - inverse of standard error
prey_f2 <- funnel(ma_prey, yaxis = "seinv", 
      xlab = "Standarized residuals",
      ylab = "Precision (inverse of SE)",  col = c(alpha("orange", 0.5)))

Figure 23— Effect size and its standard error

Figure 24— Effect size and its inverse standard error

Figure 25— Funnel plot of effect size and its standard error

Egger’s test

Code

prey_df_bias <- dt_prey %>% mutate(sqrt_inv_e_n = sqrt((Cn + Tn)/(Cn * Tn)))

prey_bias_model <- rma.mv(yi = lnRR,
                      V = lnRR_var, 
                      mods = ~ 1 + sqrt_inv_e_n,
                      random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = prey_df_bias)

summary(prey_bias_model)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-52.5371  105.0741  117.0741  134.8930  117.6873   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0921  0.3034     15     no           Study_ID 
sigma^2.2  0.0203  0.1424     58     no  Shared_control_ID 
sigma^2.3  0.0057  0.0756    123     no          Cohort_ID 
sigma^2.4  0.0323  0.1796    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 144) = 894.8647, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 144) = 0.2386, p-val = 0.6259

Model Results:

              estimate      se    tval   df    pval    ci.lb   ci.ub    
intrcpt         0.1831  0.1390  1.3172  144  0.1899  -0.0917  0.4579    
sqrt_inv_e_n    0.2307  0.4722  0.4885  144  0.6259  -0.7026  1.1640    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(prey_bias_model,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size")

Figure 26— Egger’s test of publication bias

Decline effect

Code

prey_year_model <- rma.mv(yi = lnRR,
                      V = VCV_prey, 
                      mods = ~  1 + Year,
                      random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = prey_df_bias)

summary(prey_year_model)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-51.4600  102.9199  114.9199  132.7388  115.5331   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0663  0.2575     15     no           Study_ID 
sigma^2.2  0.0225  0.1500     58     no  Shared_control_ID 
sigma^2.3  0.0045  0.0673    123     no          Cohort_ID 
sigma^2.4  0.0330  0.1818    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 144) = 729.7948, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 144) = 2.4395, p-val = 0.1205

Model Results:

         estimate       se     tval   df    pval     ci.lb     ci.ub    
intrcpt   54.3365  34.6420   1.5685  144  0.1190  -14.1362  122.8091    
Year      -0.0269   0.0172  -1.5619  144  0.1205   -0.0610    0.0071    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(prey_year_model,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication")

Figure 27— Decline effect

multi-moderator model

Code

prey_bias_all <- rma.mv(yi = lnRR,
                      V = VCV_prey, 
                      mods = ~ 1 + sqrt_inv_e_n + Year,
                      random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = prey_df_bias)

summary(prey_bias_all)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-50.9612  101.9225  115.9225  136.6624  116.7521   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0703  0.2651     15     no           Study_ID 
sigma^2.2  0.0229  0.1513     58     no  Shared_control_ID 
sigma^2.3  0.0046  0.0679    123     no          Cohort_ID 
sigma^2.4  0.0330  0.1816    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 143) = 728.4019, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 143) = 1.1593, p-val = 0.3166

Model Results:

              estimate       se     tval   df    pval     ci.lb     ci.ub    
intrcpt        53.7172  36.2602   1.4814  143  0.1407  -17.9581  125.3925    
sqrt_inv_e_n    0.0224   0.4725   0.0475  143  0.9622   -0.9116    0.9565    
Year           -0.0266   0.0180  -1.4762  143  0.1421   -0.0622    0.0090    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.5 Additional analyses

3.5.1 Background area (mm²)

Code

prey_mr_background <- rma.mv(yi = lnRR,
                        V = VCV_prey, 
                        mods = ~ Log_background,
                        random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                        test = "t",
                        method = "REML",
                        sparse = TRUE,
                        data = dt_prey)

summary(prey_mr_background)

Multivariate Meta-Analysis Model (k = 146; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-52.0900  104.1801  116.1801  133.9990  116.7932   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0825  0.2872     15     no           Study_ID 
sigma^2.2  0.0208  0.1443     58     no  Shared_control_ID 
sigma^2.3  0.0052  0.0724    123     no          Cohort_ID 
sigma^2.4  0.0329  0.1814    146     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 144) = 897.2562, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 144) = 1.0448, p-val = 0.3084

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt          -0.4706  0.6938  -0.6782  144  0.4987  -1.8420  0.9009    
Log_background    0.1048  0.1026   1.0222  144  0.3084  -0.0979  0.3075    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(prey_mr_background,
            mod = "Log_background",
            group = "Study_ID",
            xlab = "Log-transformed backgroud")

Figure 28— Effect size and background area

3.6 Figure galary

:::{.panel-tabset}

3.6.1 Discrete moderators

Code

# overall effect
prey_p1 <- orchard_plot(ma_prey,
                    group = "Study_ID",
                    xlab = "log response ratio (lnRR)", 
                    angle = 45,
                    legend.pos = "bottom.right") + 
                    scale_x_discrete(labels = c("Overall effect")) + 
                    theme(legend.direction = "horizontal",
                    legend.text = element_text(size = 6)) +
                    
                    scale_colour_okabe_ito()+
                    scale_fill_okabe_ito()

#eyespot vs. conspicuous patterns
prey_p2 <- orchard_plot(prey_mr_eyespot,
                  mod = "Treatment_stimulus",
                  group = "Study_ID",
                  xlab = "log response ratio (lnRR)",
                  angle = 45,
                  legend.pos = "bottom.right") +
                  theme(legend.direction = "horizontal",
                  legend.text = element_text(size = 6)) +
                  scale_colour_okabe_ito(order = 2:3)+
                  scale_fill_okabe_ito(order = 2:3)
                  
# prey type - real vs. artificial
prey_p3 <- orchard_plot(prey_mr_prey_type,
                  mod = "Type_prey",
                  group = "Study_ID",
                  xlab = "log response ratio (lnRR)",
                  angle = 45,
                  legend.pos = "bottom.right") +
                  theme(legend.direction = "horizontal",
                  legend.text = element_text(size = 6))+ 
                  scale_colour_okabe_ito(order = 4:5) +
                  scale_fill_okabe_ito(order = 4:5)

# combine
prey_patch1 <- prey_p1 | (prey_p2 / prey_p3)
prey_patch1 + plot_annotation(tag_levels = "a")

# output figure as pdf
# ggsave("fig1.pdf", width = 10, height = 5, dpi = 450)

Figure 29— Discrete moderators

3.6.2 Continuous moderators

Code

# these figs were created by multi meta-regression model results
# log-transformed area
prey_p4 <- bubble_plot(prey_mr_area,
                  mod = "Log_area",
                  group = "Study_ID",
                  xlab = "Log-transformed area") +
                  scale_x_continuous(breaks = seq(0,7,1.5)) +
                  scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5))

# number of patterns
prey_p5 <- bubble_plot(prey_mr_num,
            mod = "Number_pattern",
            group = "Study_ID",
            xlab = "Number of patterns") +
                  scale_x_continuous(breaks = seq(0,11,1.5)) +
                  scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5))

# combine
prey_p4 / prey_p5 + plot_annotation(tag_levels = "a")

# output figure as pdf
# ggsave("fig2_multi.pdf", width = 10, height = 10, dpi = 450)

Figure 30— Continuous moderators

3.6.3 Publication bias

Code

# funnel plot

# pdf("fig3.pdf")

funnel(ma_prey, yaxis = "seinv", 
      xlab = "Standarised residuals",
      ylab = "Precision (inverse of SE)",
      xlim = c(-4.0, 4.5), 
      ylim = c(0.01, 80.0), 
      col = c(alpha("mediumvioletred", 0.5)),
      cex = 0.7)

# dev.off()

# Egger's test and decline effect
prey_p7 <- bubble_plot(prey_bias_model,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size") +
            scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5)) 

prey_p8 <- bubble_plot(prey_year_model,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication") +  scale_y_continuous(breaks = seq(-2.5, 4.0, 1.5)) 

# combine plots created by bubble_plot()
prey_pub <- prey_p7 / prey_p8
prey_pub +  plot_annotation(tag_levels = 'a') 

# output figure as pdf
# ggsave("fig4.pdf", width = 10, height = 10, dpi = 450)

(a) Funnel plot

(b) Egger’s test and Decline effect

Figure 31— Publication bias

4 Meta-analysis - all dataset

We used meta-analytical models to calculate total I2 (a measure of heterogeneity not caused by sampling error; Higgins. et al., 2003) and the partial I2 explained by each random factor.

Which random effects remove in our models?

We can remove Shared control ID and Cohort ID

Code

# I exclude cohort_ID because sigma^2.2 = 0 and I2 = 0
ma_all <- rma.mv(yi = lnRR,
                  V = VCV, 
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_all)

summary(ma_all)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-248.5672   497.1344   507.1344   524.9761   507.3688   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0495  0.2224     32     no           Study_ID 
sigma^2.2  0.0000  0.0000     88     no  Shared_control_ID 
sigma^2.3  0.0000  0.0000    157     no          Cohort_ID 
sigma^2.4  0.2433  0.4932    263     no             Obs_ID 

Test for Heterogeneity:
Q(df = 262) = 2726.4404, p-val < .0001

Model Results:

estimate      se    tval   df    pval   ci.lb   ci.ub      
  0.2077  0.0600  3.4593  262  0.0006  0.0895  0.3259  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

i2 <- round(i2_ml(ma_all), 4)
r2 <- round(r2_ml(ma_all), 4)

I2_Total	I2_Study_ID	I2_Shared_control_ID	I2_Cohort_ID	I2_Obs_ID
96.6694	16.3369	0	0	80.3326

R2_marginal	R2_conditional
0	0.169

Table 3 - Heterogeneity of effect size

Code

orchard_plot(ma_all,
              group = "Study_ID",
              xlab = "log response ratio (lnRR)", 
              angle = 45) + 
              scale_x_discrete(labels = c("Overall effect")) + 
              scale_colour_okabe_ito()+
              scale_fill_okabe_ito()

Figure 32— Estimates of overall effect size and 95% confidence intervals

4.1 Meta-regressions: uni-moderator

Pattern type

Eyespot vs. conspicuous patterns - Is there significant difference of effect size between two patterns?

Code

#normal model
mr_eyespot <- rma.mv(yi = lnRR,
                    V = VCV, 
                    mods = ~ Treatment_stimulus,
                    random = list(~1 | Study_ID,
                                  ~1 | Obs_ID),
                    test = "t",
                    method = "REML", 
                    sparse = TRUE,
                    data = dt_all)

summary(mr_eyespot)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-247.7885   495.5769   503.5769   517.8350   503.7332   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0539  0.2321     32     no  Study_ID 
sigma^2.2  0.2436  0.4935    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 261) = 2715.7394, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 261) = 0.1562, p-val = 0.6930

Model Results:

                           estimate      se    tval   df    pval    ci.lb 
intrcpt                      0.1779  0.0984  1.8082  261  0.0717  -0.0158 
Treatment_stimuluseyespot    0.0435  0.1100  0.3953  261  0.6930  -0.1732 
                            ci.ub    
intrcpt                    0.3717  . 
Treatment_stimuluseyespot  0.2601    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_1 <- round(r2_ml(mr_eyespot), 4)

R2_marginal	R2_conditional
0.0013	0.1822

Table XX. Model goodness-of-fit

Code

orchard_plot(mr_eyespot,
            mod = "Treatment_stimulus",
            group = "Study_ID",
            xlab = "log response ratio (lnRR)",
            angle = 45) + 
            scale_colour_okabe_ito(order = 2:3)+
            scale_fill_okabe_ito(order = 2:3)

Figure 33— Eyespot vs. conspicuous patterns

Code

# intercept-removed model
mr_eyespot1 <- rma.mv(yi = lnRR,
                      V = VCV, 
                      mods = ~ Treatment_stimulus -1,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = dt_all)

summary(mr_eyespot1)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-247.7885   495.5769   503.5769   517.8350   503.7332   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0539  0.2321     32     no  Study_ID 
sigma^2.2  0.2436  0.4935    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 261) = 2715.7394, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 261) = 5.8135, p-val = 0.0034

Model Results:

                               estimate      se    tval   df    pval    ci.lb 
Treatment_stimulusconspicuous    0.1779  0.0984  1.8082  261  0.0717  -0.0158 
Treatment_stimuluseyespot        0.2214  0.0699  3.1679  261  0.0017   0.0838 
                                ci.ub     
Treatment_stimulusconspicuous  0.3717   . 
Treatment_stimuluseyespot      0.3590  ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pattern diameter (mm)

Code

mr_diameter <- rma.mv(yi = lnRR,
                      V = VCV, 
                      mods = ~ Log_diameter,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = dt_all)

summary(mr_diameter)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-244.7875   489.5749   497.5749   511.8330   497.7312   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0448  0.2116     32     no  Study_ID 
sigma^2.2  0.2383  0.4881    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 261) = 2716.5292, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 261) = 5.9052, p-val = 0.0158

Model Results:

              estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt        -0.1667  0.1643  -1.0143  261  0.3114  -0.4903  0.1569    
Log_diameter    0.1937  0.0797   2.4301  261  0.0158   0.0367  0.3507  * 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_2 <- round(r2_ml(mr_diameter), 4)

R2_marginal	R2_conditional
0.0658	0.2137

Table XX. Model goodness-of-fit

Code

bubble_plot(mr_diameter,
            mod = "Log_diameter",
            group = "Study_ID",
            xlab = "Log-transformed diameter")

Pattern area (mm²)

Code

mr_area <- rma.mv(yi = lnRR,
                  V = VCV, 
                  mods = ~ Log_area,
                  random = list(~1 | Study_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_all)

summary(mr_area)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-244.2677   488.5354   496.5354   510.7935   496.6916   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0491  0.2216     32     no  Study_ID 
sigma^2.2  0.2362  0.4860    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 261) = 2713.1669, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 261) = 6.8516, p-val = 0.0094

Model Results:

          estimate      se     tval   df    pval    ci.lb   ci.ub     
intrcpt    -0.1841  0.1610  -1.1437  261  0.2538  -0.5010  0.1329     
Log_area    0.1103  0.0421   2.6176  261  0.0094   0.0273  0.1933  ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_3 <- round(r2_ml(mr_area), 4)

R2_marginal	R2_conditional
0.0816	0.2396

Table XX. Model goodness-of-fit

Code

bubble_plot(mr_area,
            mod = "Log_area",
            group = "Study_ID",
            xlab = "Log-transformed area")

Pattern number

Code

mr_num <- rma.mv(yi = lnRR,
                  V = VCV,
                  mods = ~ Number_pattern,
                  random = list(~1 | Study_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_all)

summary(mr_num)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-245.9910   491.9821   499.9821   514.2402   500.1383   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0473  0.2176     32     no  Study_ID 
sigma^2.2  0.2395  0.4894    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 261) = 2703.3403, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 261) = 4.4221, p-val = 0.0364

Model Results:

                estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt           0.3444  0.0881   3.9111  261  0.0001   0.1710   0.5178  *** 
Number_pattern   -0.0578  0.0275  -2.1029  261  0.0364  -0.1120  -0.0037    * 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_4 <- round(r2_ml(mr_num), 4)

R2_marginal	R2_conditional
0.0195	0.1813

Table XX. Model goodness-of-fit

Code

bubble_plot(mr_num,
            mod = "Number_pattern",
            group = "Study_ID",
            xlab = "Number of patterns") 

Prey material type

Our dataset includes two material types of prey: live and artificial. Is there significant difference of effect size between two stimuli?

Code

# normal model
mr_prey_type <- rma.mv(yi = lnRR,
                  V = VCV, 
                  mods = ~ Type_prey,
                  random = list(~1 | Study_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_all)

summary(mr_prey_type)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-247.6124   495.2249   503.2249   517.4830   503.3811   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0571  0.2390     32     no  Study_ID 
sigma^2.2  0.2424  0.4923    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 261) = 2687.4935, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 261) = 0.3362, p-val = 0.5625

Model Results:

               estimate      se    tval   df    pval    ci.lb   ci.ub    
intrcpt          0.1834  0.0762  2.4061  261  0.0168   0.0333  0.3334  * 
Type_preyreal    0.0772  0.1331  0.5798  261  0.5625  -0.1849  0.3393    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_5 <- round(r2_ml(mr_prey_type), 4)

R2_marginal	R2_conditional
0.0035	0.1935

Table XX. Model goodness-of-fit

Code

orchard_plot(mr_prey_type,
              mod = "Type_prey",
              group = "Study_ID",
              xlab = "Type of prey",
              angle = 45) +
            scale_colour_okabe_ito(order = 4:5)+
            scale_fill_okabe_ito(order = 4:5)

Code

# intercept-removed model
mr_prey_type1 <- rma.mv(yi = lnRR,
                        V = VCV, 
                        mods = ~ Type_prey -1,
                        random = list(~1 | Study_ID,
                                      ~1 | Obs_ID),
                        test = "t",
                        method = "REML", 
                        sparse = TRUE,
                        data = dt_all)

summary(mr_prey_type1)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-247.6124   495.2249   503.2249   517.4830   503.3811   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0571  0.2390     32     no  Study_ID 
sigma^2.2  0.2424  0.4923    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 261) = 2687.4935, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 261) = 5.7450, p-val = 0.0036

Model Results:

                     estimate      se    tval   df    pval   ci.lb   ci.ub    
Type_preyartificial    0.1834  0.0762  2.4061  261  0.0168  0.0333  0.3334  * 
Type_preyreal          0.2605  0.1091  2.3876  261  0.0177  0.0457  0.4754  * 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prey shape type

Our dataset includes 4 types of prey type: real butterfly, abstract butterfly, abstract caterpillar, and abstract prey. Is there significant difference of effect size between two stimuli?

Code

mr_prey_shape <- rma.mv(yi = lnRR,
                        V = VCV, 
                        mods = ~ Shape_prey -1,
                        random = list(~1 | Study_ID,
                                      ~1 | Obs_ID),
                                      test = "t",
                                      method = "REML", 
                                      sparse = TRUE,
                                      data = dt_all)

summary(mr_prey_shape)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-243.8328   487.6657   499.6657   521.0067   499.9990   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0562  0.2371     32     no  Study_ID 
sigma^2.2  0.2399  0.4898    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 259) = 2494.5220, p-val < .0001

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 259) = 3.7295, p-val = 0.0057

Model Results:

                                estimate      se    tval   df    pval    ci.lb 
Shape_preyabstract_butterfly      0.3223  0.1080  2.9841  259  0.0031   0.1096 
Shape_preyabstract_caterpillar    0.0663  0.1286  0.5160  259  0.6063  -0.1868 
Shape_preyabstract_prey           0.0126  0.1879  0.0671  259  0.9465  -0.3574 
Shape_preybutterfly               0.2600  0.1085  2.3964  259  0.0173   0.0463 
                                 ci.ub     
Shape_preyabstract_butterfly    0.5349  ** 
Shape_preyabstract_caterpillar  0.3195     
Shape_preyabstract_prey         0.3826     
Shape_preybutterfly             0.4737   * 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_6 <- round(r2_ml(mr_prey_shape), 4)

dat_all$Shape_prey <- factor(dat_all$Shape_prey)
levels(dat_all$Shape_prey)

[1] "abstract_butterfly"   "abstract_caterpillar" "abstract_prey"       
[4] "butterfly"           

Code

mat_ex <- 
    cbind(contrMat(rep(1,
                       length(mr_prey_shape$ci.ub)),
                   type = "Tukey"))
sig_test <- summary(glht(mr_prey_shape, linfct= mat_ex), test = adjusted("none"))

sig_test

     Simultaneous Tests for General Linear Hypotheses

Fit: rma.mv(yi = lnRR, V = VCV, mods = ~Shape_prey - 1, random = list(~1 | 
    Study_ID, ~1 | Obs_ID), data = dt_all, method = "REML", test = "t", 
    sparse = TRUE)

Linear Hypotheses:
           Estimate Std. Error z value Pr(>|z|)
2 - 1 == 0 -0.25592    0.16791  -1.524    0.127
3 - 1 == 0 -0.30966    0.21673  -1.429    0.153
4 - 1 == 0 -0.06227    0.15309  -0.407    0.684
3 - 2 == 0 -0.05373    0.22768  -0.236    0.813
4 - 2 == 0  0.19365    0.16823   1.151    0.250
4 - 3 == 0  0.24739    0.21698   1.140    0.254
(Adjusted p values reported -- none method)

R2_marginal	R2_conditional
0.0542	0.2338

Table XX. Model goodness-of-fit

Code

orchard_plot(mr_prey_shape,
              mod = "Shape_prey",
              group = "Study_ID",
              xlab = "Shape of prey",
              angle = 45) +
              scale_colour_okabe_ito(order = 6:9)+
              scale_fill_okabe_ito(order = 6:9)

Figure 34— Effect size and prey shape

4.2 Correlation visualisation and choose moderators

Before we run multi-moderator meta-regressions, we need to consider the correlation between moderators. Area, diameter and background seem to be correlated. We visualised the correlation between these variables.

Visualisation of correlation

Code

# correlation between continuous variables
corr_cont <- round(cor(dt_all[, c("Diameter_pattern", "Area_pattern", "Number_pattern", "Area_background")]),2)

p1 <- ggcorrplot(corr_cont, hc.order = TRUE, lab = TRUE, outline.col = "white", colors = c("#6D9EC1", "white", "#E46726"), title = "(a) Continuous variables")

corr_cont_log <- round(cor(dt_all[, c("Log_diameter", "Log_area", "Number_pattern", "Log_background")]),2)

p2 <- ggcorrplot(corr_cont_log, hc.order = TRUE, lab = TRUE, outline.col = "white", colors = c("#6D9EC1", "white", "#E46726"), title = "(b) Log-transormed continuous variables")

p1 + p2 +plot_layout(guides = 'collect')

Figure 35— Correlation between coninuous moderators

Code

dat1 <- dt_all %>%
  dplyr::select("Treatment_stimulus", "Type_prey", "Shape_prey")

corr_cat <- GKtauDataframe(dat1)
plot(corr_cat)

Figure 36— Correlation between categorical variables

We should not include “Shape_prey” and “Type_prey” in the model at the same time. Because they are correlated.

Choose moderators

We used model R2 values to find better model and modelator VIF values to check multicollinearity between moderators. Higher R2 indicates better model and VIF > 2 indicates multicollinearity.

4.2.0.1 R2

Code

r2_area <- rma.mv(yi = lnRR,
                  V = VCV,
                  mods = ~ Log_area + Number_pattern,
                  random = list(~1 | Study_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML",
                  sparse = TRUE,
                  data = dt_all)

r2_diameter <- rma.mv(yi = lnRR,
                      V = VCV,
                      mods = ~ Log_diameter + Number_pattern,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_all)

# check r2 values
r2_ml(r2_area)

   R2_marginal R2_conditional 
    0.07553592     0.22788272 

Code

r2_ml(r2_diameter)

   R2_marginal R2_conditional 
    0.06439709     0.20868706 

It seems area is a better predictor than diameter .

4.3 Meta-regressions: multi-moderator

All potential moderators (full model)

Include all moderators in the model.

Code

mr_all <- rma.mv(yi = lnRR,
                V = VCV, 
                mods = ~ Treatment_stimulus + Log_area 
                + Number_pattern + Type_prey,
                random = list(~1 | Study_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_all)

summary(mr_all)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-240.5439   481.0877   495.0877   519.9584   495.5357   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0555  0.2355     32     no  Study_ID 
sigma^2.2  0.2342  0.4840    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 258) = 2592.2639, p-val < .0001

Test of Moderators (coefficients 2:5):
F(df1 = 4, df2 = 258) = 2.6696, p-val = 0.0327

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.0883  0.2151  -0.4104  258  0.6818  -0.5117 
Treatment_stimuluseyespot    0.0103  0.1144   0.0897  258  0.9286  -0.2150 
Log_area                     0.0979  0.0457   2.1410  258  0.0332   0.0079 
Number_pattern              -0.0504  0.0300  -1.6802  258  0.0941  -0.1095 
Type_preyreal                0.1924  0.1438   1.3379  258  0.1821  -0.0908 
                            ci.ub    
intrcpt                    0.3352    
Treatment_stimuluseyespot  0.2355    
Log_area                   0.1879  * 
Number_pattern             0.0087  . 
Type_preyreal              0.4756    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(mr_all)

   R2_marginal R2_conditional 
    0.07939435     0.25570867 

Moderators associated with conspicuousness

Only include moderators related to conspicuousness (pattern area and number of patterns) in the model.

Code

mr_cons <- rma.mv(yi = lnRR,
                V = VCV, 
                mods = ~ Log_area + Number_pattern,
                random = list(~1 | Study_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_all)

summary(mr_cons)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.9378   485.8757   495.8757   513.6791   496.1119   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0466  0.2158     32     no  Study_ID 
sigma^2.2  0.2360  0.4858    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 260) = 2679.1464, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 260) = 4.3867, p-val = 0.0134

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt          -0.0218  0.1969  -0.1109  260  0.9118  -0.4095  0.3659    
Log_area          0.0908  0.0438   2.0735  260  0.0391   0.0046  0.1770  * 
Number_pattern   -0.0394  0.0287  -1.3742  260  0.1706  -0.0960  0.0171    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(mr_cons)

   R2_marginal R2_conditional 
    0.07553592     0.22788272 

Code

bubble_plot(mr_cons,
            mod = "Log_area",
            group = "Study_ID",
            xlab = "Log-transformed area")

Code

mr_cons_1 <- rma.mv(yi = lnRR,
                V = VCV, 
                mods = ~ Log_area + Number_pattern + Treatment_stimulus,
                random = list(~1 | Study_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_all)

summary(mr_cons_1)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.1568   484.3135   496.3135   517.6545   496.6469   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0513  0.2264     32     no  Study_ID 
sigma^2.2  0.2361  0.4859    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 259) = 2678.9397, p-val < .0001

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 259) = 2.9683, p-val = 0.0325

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.0563  0.2109  -0.2668  259  0.7899  -0.4715 
Log_area                     0.0914  0.0448   2.0403  259  0.0423   0.0032 
Number_pattern              -0.0404  0.0290  -1.3949  259  0.1642  -0.0975 
Treatment_stimuluseyespot    0.0517  0.1085   0.4764  259  0.6342  -0.1620 
                            ci.ub    
intrcpt                    0.3590    
Log_area                   0.1796  * 
Number_pattern             0.0166    
Treatment_stimuluseyespot  0.2653    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(mr_cons_1)

   R2_marginal R2_conditional 
    0.07948479     0.24369201 

Moderators associated with prey

Code

mr_pre <- rma.mv(yi = lnRR,
                V = VCV, 
                mods = ~ Treatment_stimulus + Type_prey,
                random = list(~1 | Study_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_all)

summary(mr_pre)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-246.8556   493.7111   503.7111   521.5145   503.9474   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0608  0.2466     32     no  Study_ID 
sigma^2.2  0.2431  0.4930    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 260) = 2667.6509, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 260) = 0.2089, p-val = 0.8116

Model Results:

                           estimate      se    tval   df    pval    ci.lb 
intrcpt                      0.1650  0.1037  1.5906  260  0.1129  -0.0393 
Treatment_stimuluseyespot    0.0300  0.1171  0.2562  260  0.7980  -0.2005 
Type_preyreal                0.0699  0.1412  0.4947  260  0.6212  -0.2082 
                            ci.ub    
intrcpt                    0.3693    
Treatment_stimuluseyespot  0.2605    
Type_preyreal              0.3480    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(mr_pre)

   R2_marginal R2_conditional 
    0.00415544     0.20344953 

4.4 Publication bias

Funnel plot

Code

# funnel plot - standard error
funnel(ma_all, yaxis = "sei", 
      xlab = "Standarised residuals",
      ylab = "Precision (inverse of SE)" )

# funnel plot - inverse of standard error
f2 <- funnel(ma_all, yaxis = "seinv", 
      xlab = "Standarized residuals",
      ylab = "Precision (inverse of SE)",  col = c(alpha("orange", 0.5)))

Figure 37— Effect size and its standard error

Figure 38— Effect size and its inverse standard error

Figure 39— Funnel plot of effect size and its standard error

Egger’s test

Code

df_bias <- dt_all %>% mutate(sqrt_inv_e_n = sqrt((Cn + Tn)/(Cn * Tn)))

bias_model <- rma.mv(yi = lnRR,
                      V = lnRR_var, 
                      mods = ~ 1 + sqrt_inv_e_n,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = df_bias)

summary(bias_model)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-246.1777   492.3554   500.3554   514.6135   500.5116   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0960  0.3099     32     no  Study_ID 
sigma^2.2  0.2168  0.4656    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 261) = 2152.6798, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 261) = 0.0084, p-val = 0.9269

Model Results:

              estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt         0.2426  0.1382   1.7550  261  0.0804  -0.0296  0.5148  . 
sqrt_inv_e_n   -0.0357  0.3884  -0.0918  261  0.9269  -0.8004  0.7291    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(bias_model,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size")

Figure 40— Egger’s test of publication bias

Decline effect

Code

year_model <- rma.mv(yi = lnRR,
                      V = VCV, 
                      mods = ~  1 + Year,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = df_bias)

summary(year_model)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-247.6851   495.3701   503.3701   517.6282   503.5264   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0525  0.2291     32     no  Study_ID 
sigma^2.2  0.2437  0.4937    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 261) = 2653.9609, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 261) = 0.0059, p-val = 0.9389

Model Results:

         estimate       se     tval   df    pval     ci.lb    ci.ub    
intrcpt    1.2075  13.0311   0.0927  261  0.9262  -24.4519  26.8670    
Year      -0.0005   0.0065  -0.0767  261  0.9389   -0.0133   0.0123    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(year_model,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication")

Figure 41— Decline effect

Multi-moderator

Code

multi_bias <- rma.mv(yi = lnRR,
                      V = VCV, 
                      mods = ~sqrt_inv_e_n + Year,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = df_bias)

summary(multi_bias)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-246.4927   492.9854   502.9854   520.7888   503.2217   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0458  0.2140     32     no  Study_ID 
sigma^2.2  0.2456  0.4956    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 260) = 2597.3634, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 260) = 0.3597, p-val = 0.6982

Model Results:

              estimate       se     tval   df    pval     ci.lb    ci.ub    
intrcpt         4.6120  13.2052   0.3493  260  0.7272  -21.3907  30.6147    
sqrt_inv_e_n   -0.3289   0.3916  -0.8400  260  0.4017   -1.1000   0.4421    
Year           -0.0021   0.0066  -0.3270  260  0.7439   -0.0151   0.0108    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(multi_bias,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size")

Code

bubble_plot(multi_bias,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication")

4.5 Additional analyses

Dataset

We compared effect size between two datasets (predator and prey) to see whether there is any difference in effect size between two datasets.

Code

mr_dataset <- rma.mv(yi = lnRR,
                      V = VCV,
                      mods = ~ Dataset,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_all)

summary(mr_dataset)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-247.3783   494.7565   502.7565   517.0146   502.9128   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0525  0.2291     32     no  Study_ID 
sigma^2.2  0.2429  0.4928    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 261) = 2649.6615, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 261) = 0.6243, p-val = 0.4302

Model Results:

             estimate      se    tval   df    pval    ci.lb   ci.ub    
intrcpt        0.1575  0.0885  1.7798  261  0.0763  -0.0167  0.3316  . 
Datasetprey    0.0944  0.1195  0.7901  261  0.4302  -0.1409  0.3297    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

orchard_plot(mr_dataset,
            mod = "Dataset",
            group = "Study_ID",
            xlab = "Dataset", 
            angle = 45) 

Figure 42— Effect size and dataset - predator and prey

Background area (mm²)

Code

mr_background <- rma.mv(yi = lnRR,
                        V = VCV, 
                        mods = ~ Log_background,
                        random = list(~1 | Study_ID,
                                      ~1 | Obs_ID),
                        test = "t",
                        method = "REML",
                        sparse = TRUE,
                        data = dt_all)

summary(mr_background)

Multivariate Meta-Analysis Model (k = 263; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-247.3566   494.7131   502.7131   516.9712   502.8694   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0559  0.2364     32     no  Study_ID 
sigma^2.2  0.2429  0.4929    263     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 261) = 2660.5568, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 261) = 0.2498, p-val = 0.6176

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt           0.4218  0.4312   0.9782  261  0.3289  -0.4273  1.2709    
Log_background   -0.0306  0.0613  -0.4998  261  0.6176  -0.1513  0.0900    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(mr_background)

   R2_marginal R2_conditional 
   0.004611716    0.190784945 

Code

bubble_plot(mr_background,
            mod = "Log_background",
            group = "Study_ID",
            xlab = "Log-transformed backgroud")

Figure 43— Effect size and background area

4.6 Figure galary

Discrete moderators

Code

# overall effect
p1 <- orchard_plot(ma_all,
                    group = "Study_ID",
                    xlab = "log response ratio (lnRR)", 
                    angle = 45,
                    legend.pos = "bottom.right") + 
                    scale_x_discrete(labels = c("Overall effect")) + 
                    theme(legend.direction = "horizontal",
                    legend.text = element_text(size = 6)) +
                    
                    scale_colour_okabe_ito()+
                    scale_fill_okabe_ito()

#eyespot vs. conspicuous patterns
p2 <- orchard_plot(mr_eyespot,
                  mod = "Treatment_stimulus",
                  group = "Study_ID",
                  xlab = "log response ratio (lnRR)",
                  angle = 45,
                  legend.pos = "bottom.right") +
                  theme(legend.direction = "horizontal",
                  legend.text = element_text(size = 6)) +
                  scale_colour_okabe_ito(order = 2:3)+
                  scale_fill_okabe_ito(order = 2:3)
                  
# prey type - real vs. artificial
p3 <- orchard_plot(mr_prey_type,
                  mod = "Type_prey",
                  group = "Study_ID",
                  xlab = "log response ratio (lnRR)",
                  angle = 45,
                  legend.pos = "bottom.right") +
                  theme(legend.direction = "horizontal",
                  legend.text = element_text(size = 6))+ 
                  scale_colour_okabe_ito(order = 4:5) +
                  scale_fill_okabe_ito(order = 4:5)

# combine
patch1 <- p1 | (p2 / p3)
patch1 + plot_annotation(tag_levels = "a")

# output figure as pdf
# ggsave("fig1.pdf", width = 10, height = 5, dpi = 450)

Figure 44— Discrete moderators

Continuous moderators

Code

# these figs were created by multi meta-regression model results
# log-transformed area
p4 <- bubble_plot(mr_cons,
                  mod = "Log_area",
                  group = "Study_ID",
                  xlab = "Log-transformed area") +
                  scale_x_continuous(breaks = seq(0,7,1.5)) +
                  scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5))

# number of patterns
p5 <- bubble_plot(mr_cons,
            mod = "Number_pattern",
            group = "Study_ID",
            xlab = "Number of patterns") +
                  scale_x_continuous(breaks = seq(0,11,1.5)) +
                  scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5))

# combine
p4 / p5 + plot_annotation(tag_levels = "a")

# output figure as pdf
# ggsave("fig2_multi.pdf", width = 10, height = 10, dpi = 450)

Figure 45— Continuous moderators

Publication bias

Code

# funnel plot

# pdf("fig3.pdf")

funnel(ma_all, yaxis = "seinv", 
      xlab = "Standarised residuals",
      ylab = "Precision (inverse of SE)",
      xlim = c(-4.0, 4.5), 
      ylim = c(0.01, 80.0), 
      col = c(alpha("mediumvioletred", 0.5)),
      cex = 0.7)

# dev.off()

# Egger's test and decline effect
p7 <- bubble_plot(bias_model,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size") +
            scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5)) 

p8 <- bubble_plot(year_model,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication") +  scale_y_continuous(breaks = seq(-2.5, 4.0, 1.5)) 

# combine plots created by bubble_plot()
pub <- p7 / p8
pub +  plot_annotation(tag_levels = 'a') 

# output figure as pdf
# ggsave("fig4.pdf", width = 10, height = 10, dpi = 450)

Figure 46— Funnel plot

Figure 47— Egger’s test and Decline effect

Figure 48— Publication bias

5 R session infomation

R version 4.3.1 (2023-06-16)

Platform: aarch64-apple-darwin20 (64-bit)

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: parallel, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: orchaRd(v.2.0), multcomp(v.1.4-25), TH.data(v.1.1-2), MASS(v.7.3-60), survival(v.3.5-5), mvtnorm(v.1.2-3), miWQS(v.0.4.4), metafor(v.4.2-0), numDeriv(v.2016.8-1.1), metadat(v.1.2-0), Matrix(v.1.6-1), patchwork(v.1.1.3), lubridate(v.1.9.2), forcats(v.1.0.0), stringr(v.1.5.0), dplyr(v.1.1.3), purrr(v.1.0.2), readr(v.2.1.4), tidyr(v.1.3.0), tibble(v.3.2.1), tidyverse(v.2.0.0), knitr(v.1.43), here(v.1.0.1), ggcorrplot(v.0.1.4), ggplot2(v.3.4.3), ggokabeito(v.0.1.0), dtplyr(v.1.3.1), DT(v.0.29), GoodmanKruskal(v.0.0.3) and ape(v.5.7-1)

loaded via a namespace (and not attached): gridExtra(v.2.3), glm2(v.1.2.1), sandwich(v.3.0-2), rlang(v.1.1.1), magrittr(v.2.0.3), compiler(v.4.3.1), mgcv(v.1.8-42), reshape2(v.1.4.4), vctrs(v.0.6.3), quantreg(v.5.97), pkgconfig(v.2.0.3), fastmap(v.1.1.1), ellipsis(v.0.3.2), backports(v.1.4.1), mcmc(v.0.9-7), labeling(v.0.4.3), pander(v.0.6.5), utf8(v.1.2.3), rmarkdown(v.2.24), tzdb(v.0.4.0), ggbeeswarm(v.0.7.2), MatrixModels(v.0.5-2), xfun(v.0.40), cachem(v.1.0.8), jsonlite(v.1.8.7), gmm(v.1.8), cluster(v.2.1.4), R6(v.2.5.1), bslib(v.0.5.1), stringi(v.1.7.12), Rsolnp(v.1.16), rlist(v.0.4.6.2), rpart(v.4.1.19), estimability(v.1.4.1), jquerylib(v.0.1.4), Rcpp(v.1.0.11), zoo(v.1.8-12), tmvmixnorm(v.1.1.1), base64enc(v.0.1-3), pacman(v.0.5.1), splines(v.4.3.1), nnet(v.7.3-19), timechange(v.0.2.0), tidyselect(v.1.2.0), rstudioapi(v.0.15.0), yaml(v.2.3.7), codetools(v.0.2-19), plyr(v.1.8.8), lattice(v.0.21-8), withr(v.2.5.0), coda(v.0.19-4), evaluate(v.0.21), tmvtnorm(v.1.5), foreign(v.0.8-84), pillar(v.1.9.0), corrplot(v.0.92), checkmate(v.2.2.0), stats4(v.4.3.1), generics(v.0.1.3), rprojroot(v.2.0.3), invgamma(v.1.1), mathjaxr(v.1.6-0), truncnorm(v.1.0-9), hms(v.1.1.3), munsell(v.0.5.0), scales(v.1.2.1), xtable(v.1.8-4), glue(v.1.6.2), emmeans(v.1.8.8), Hmisc(v.5.1-0), tools(v.4.3.1), data.table(v.1.14.8), SparseM(v.1.81), matrixNormal(v.0.1.1), grid(v.4.3.1), MCMCpack(v.1.6-3), crosstalk(v.1.2.0), colorspace(v.2.1-0), nlme(v.3.1-162), beeswarm(v.0.4.0), condMVNorm(v.2020.1), vipor(v.0.4.5), htmlTable(v.2.4.1), latex2exp(v.0.9.6), Formula(v.1.2-5), cli(v.3.6.1), fansi(v.1.0.4), gtable(v.0.3.4), sass(v.0.4.7), digest(v.0.6.33), farver(v.2.1.1), htmlwidgets(v.1.6.2), htmltools(v.0.5.6) and lifecycle(v.1.0.3)